DESIGN OF PARAMETRICALLY FORCED PATTERNS AND 
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Abstract. The Faraday wave experiment is a classic example of a system driven by paramet- 
ric forcing, and it produces a wide range of complex patterns, including superlattico patterns and 
quasipatterns. Nonlinear three-wave interactions between driven and weakly damped modes play a 
key role in determining which patterns are favoured. We use this idea to design single and multi- 
frequency forcing functions that produce examples of superlattice patterns and quasipatterns in a 
new model PDE with parametric forcing. We make quantitative comparisons between the predicted 
patterns and the solutions of the PDE. Unexpectedly, the agreement is good only for parameter 
values very close to onset. The reason that the range of validity is limited is that the theory re- 
quires strong damping of all modes apart from the driven pattern-forming modes. This is in conflict 
with the requirement for weak damping if three- wave coupling is to influence pattern selection effec- 
tively. We distinguish the two different ways that three-wave interactions can be used to stabilise 
quasipatterns, and present examples of 12-, 14- and 20-fold approximate quasipatterns. We identify 
which computational domains provide the most accurate approximations to 12-fold quasipatterns, 
and systematically investigate the Fourier spectra of the most accurate approximations. 
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day waves. 
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1. Introduction. The classic Faraday wave experiment consists of a horizontal 
layer of fluid that spontaneously develops a pattern of standing waves on its sur- 
face as it is driven by vertical oscillation with amplitude exceeding a critical value; 
see [2, 32, 43] for recent reviews and surveys. Many other experimental, biological 
and environmental systems also form patterns [13,30], but Faraday wave experiments 
have consistently produced patterns with remarkably high degrees of symmetry. One 
consequence of this is that, over the years, Faraday wave experiments have repeatedly 
produced new patterns of behaviour that went beyond the existing range of theoretical 
understanding and required the development of new ideas for their explanation. For 
example, in the early 1990's, quasipatterns were discovered in two different Faraday 
wave experiments, one with a low-viscosity deep layer of fluid with single-frequency 
forcing [5-7], and the other with a high- viscosity shallow layer of fluid and forcing 
with two commensurate temporal frequencies [19]. These patterns are periodic in 
time but are quasiperiodic in any spatial direction, that is, the amplitude of the pat- 
tern (taken along any direction in the plane) can be regarded as the sum of waves with 
incommensurate spatial frequencies. In spite of this, the spatial Fourier transforms 
of quasipatterns have 8, 10 or 12-fold rotational order. Quasipatterns are of course 
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related to quasicrystals [53], and quasipatterns have been found in nonlinear optical 
systems [29], in shaken convection [50,64] and in liquid crystals [36] as well as the 
Faraday wave experiment [2,5-7,19,33]. There is, as yet, no satisfactory theoreti- 
cal understanding of the formation of quasipatterns owing to the problem of small 
divisors [51]. 

Theoretical efforts aimed at understanding the pattern selection problem have 
centered around weakly nonlinear theory. The calculations for the real Faraday wave 
problem, with finite depth and non-zero viscosity, are difficult (these involve solving 
the Navier-Stokes equations with a free surface boundary condition [57]). Most cal- 
culations aimed at producing superlattice patterns and quasipatterns have focussed 
on simpler equations, such as the Zhang- Viiials [67] equations, which model Navier- 
Stokes when the depth is infinite and the viscosity is small, or on model equations, such 
as variations on the Swift-Hohenberg equation [22.37,42] or the Fitzhugh-Nagumo 
equations [15]. 

Nonlinear three-wave resonant interactions have long been recognised as playing 
a key role in pattern selection in Faraday wave experiments, or other situations where 
complex patterns are found [19,41,45,68]. A series of papers [47-49,55,56,61,62] has 
developed this idea, using symmetry considerations to understand pattern selection 
in Faraday wave experiments with two-frequency forcing, exploiting the three-wave 
resonant interactions in the context of weakly broken Hamiltonian structure. This 
approach was able to explain several of the experimentally observed superlattice pat- 
terns, and suggested ways of designing multi-frequency forcing functions that could be 
used to control which patterns would emerge [48,61,62]. The approach is in principle 
predictive, but has only been used to determine which additional frequencies to add 
to the forcing function in order to make observed patterns more robust [2,16,20]. The 
theory has not been tested quantitatively against solutions of a pattern-forming sys- 
tem, though weakly nonlinear coefficients have been computed for the Zhang- Vihals 
equations [49,61] and, more recently, for the Navier-Stokes equations in the infinite 
depth case [57]. 

The main goal of this paper is to come to a greater understanding of some of the 
complex patterns that are found in high-precision large aspect ratio Faraday wave 
experiments. Much of the complexity arises from using two or more frequencies in 
the forcing of the experiment, and recent work [49, 61] explains in principle how to 
connect the amplitudes and phases of the various components of the forcing frequen- 
cies to the nonlinear pattern selection problem. The existing theory provides rules of 
thumb for designing forcing functions that should encourage the appearance of par- 
ticular patterns, and it seems to work well, at least qualitatively, and at least in some 
circumstances [2,16,20]. 

We take the point of view that in order to claim convincingly that we under- 
stand the pattern selection process in these problems with multi-frequency forcing, 
we should be able to predict in advance which patterns will be found for different 
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parameter values, we should be able to predict the amplitudes and range of stability 
of the patterns, and and we should test against a pattern forming system that is not 
constrained to produce only a limited range of patterns. In order to have a flexi- 
ble framework for testing predictions, we have devised a partial differential equation 
(PDE) with multi-frequency forcing (|3.ip that shares many of the characteristics of 
the real Faraday wave experiment, but that has easily controllable dissipation and 
dispersion relations, and simple nonlinear terms. The linear behaviour of the PDE 
reduces to the damped Mathieu equation, with subharmonic and harmonic tongues, 
and the simple quadratic and cubic nonlinearities allow three-wave interactions. The 
PDE has a Hamiltonian limit, but it differs from the real situation in the details of its 
dispersion relation, and the lack of any coupling to a large-scale mean flow. Notwith- 
standing these differences, the PDE allows us to explore in detail some of the generic 
issues surrounding pattern selection in very large aspect pattern forming systems with 
parametric forcing. A preliminary discussion of the question of quasipattern selection 
in the new model PDE can be found in [52]. 

In section [21 we review the details of how resonant three-wave interactions in- 
fluence pattern selection. The main idea is that two pattern-forming modes, with 
wavevectors separated by an angle 9, are coupled to a weakly damped mode, and 
this coupling can lead to the angle 6 either featuring in the resulting pattern or be- 
ing eliminated from the resulting pattern [19,41,45,68]. We introduce the model 
PDE in section [21 and describe its linear and nonlinear features in sections [H and [SI 
Appendix VK\ gives full details of the weakly nonlinear calculations. 

In section [6] we devise a forcing function that stabilises the 22° superlattice pat- 
terns that have been observed in large aspect ratio Faraday wave experiments [2,20, 
33] . We compute fully nonlinear solutions of the PDE and compare the computed pat- 
tern amplitudes with those predicted by weakly nonlinear theory. This demonstrates 
that the agreement is quantitatively accurate only very close to onset (within 0.1%). 
We show that the limit on the range of validity is because of the presence of the weakly 
damped modes that are required for the superlattice pattern to be stabilised. Thus 
parameter regimes that are likely to produce the most interesting patterns, arising 
from coupling to weakly damped modes, are also parameter regimes where weakly 
nonlinear theory has the most restricted validity. 

Two mechanisms have been proposed for quasipattern formation, both building 
on ideas of Newell and Pomeau [45], and one aim of this paper is to demonstrate 
that both proposed mechanisms for quasipattern formation are viable (preliminary 
work is reported in [52]). One mechanism applies to single frequency forced Faraday 
waves [67] and has been tested experimentally [66] . Another was developed to explain 
the origin of the two length scales in superlattice patterns [49, 62] found in two- 
frequency experiments [33]. The ideas have not been tested quantitatively, but have 
been used qualitatively to control quasipattern [2, 16] and superlattice pattern [20] 
formation in two and three- frequency experiments. We explore the two mechanisms 
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of quasipattcrn formation in the model PDE in sections [7] and [SI 

Also in section [3 we address the distinction between true and approximate quasi- 
patterns, as found in numerical experiments with periodic boundary conditions. Ow- 
ing to the problem of small divisors, there is as yet no satisfactory mathematical 
treatment of quasipatterns [51]. In spite of this, the weakly nonlinear stability cal- 
culations, which are in the framework of a 12-mode amplitude expansion truncated 
at cubic order, prove to be a reliable guide to finding parameter values where ap- 
proximate quasipatterns are stable. The fact that stable 12-fold quasipatterns are 
found where they are expected demonstrates that this approach provides useful in- 
formation, in spite of the reservations expressed in [51]. We explore the effect of 
domain size on the accuracy of the approximation to a true quasipattcrn, and show 
how certain domains yield particularly accurate approximations. 

In section [51 we present examples of turbulent crystals [45] : situations in which 
Fourier modes oriented more than about 20° apart do not affect each other, at the 
level of a cubic truncation. In this context, we find examples of 12-, 14- and 20-fold 
quasipatterns. These are the first examples of quasipatterns of order greater than 12 
found as stable solutions of a PDE (a preliminary presentation of the 14-fold example 
is in [52]). 

We summarise our result in section [9] 

2. Theoretical background. Resonant triads have played a key role in our 
understanding of pattern formation [19,41,45,68]. This section reviews, somewhat 
heuristically, the basic selection mechanisms in the case of wave patterns that are 
parametrically pumped by a two (or more) frequency forcing function. The details 
behind this summary can be found in [48,61]. We write the forcing function as 

f{t) = /m COs(mt + Cf)rn) + fn COs{nt + 0„) + (2.1) 

where m and n are integers, fm and /„ are real amplitudes, and (j)m and (/>„ are 
phases. (We could, of course, set (pm ~ without loss of generality.) Here we consider 
m to be the dominant driving frequency, and focus on a pair of waves, each with 
wavenumber fc„j, which satisfies the dispersion relation ri(fc,„) ~ m/2 associated with 
the linearized problem. In other words, these waves naturally oscillate at a frequency 
that is subharmonic to the dominant driving frequency m, and are typically the eas- 
iest to excite parametrically. We write the critical modes in the form 2:]^e*'°^'^"'"*'"*/^ 
and 2;2e*'"^ '^+""*/^ (together with their complex conjugates), neglecting the higher 
temporal frequency contributions to the waves. These waves will interact nonlinearly 
with waves ijue^i'3-x+in{k3)t ^ where u; is a complex amplitude, fcg = fci +^2 and 0(^3) 
is the frequency associated with k^, provided that either (1) the same resonance con- 
dition is met with the temporal frequencies, i.e., ^I{k3) = ^ + ^, as in figure [TlT a.b). 
or (2) any mismatch A = |r2(fc3) — ^ ~ tI ^-'^^^ temporal resonance condition can 
be compensated for by the forcing function f{t) (figure I^TTb ). Case (1) corresponds to 
1 : 2 resonance, which occurs even for single frequency forcing (/„ ~ 0), and case (2) 




Fig. 2.1. (a) 1 : 2 resonance occurs between two modes with wavevectors fci and (blue), 
with the same wavenumber fci and separated by an angle 6, and a third mode with wavevector ks 
(magenta), provided fca = k±+k2 and Q{k:j) = 2Q{ki). (b) A special case o/l : 2 resonance in space 
and time occurs for 8 = when the dispersion relation satisfies Q{2ki) = 2Q{ki). (c) With two- 
freguency m : n forcing, the nonlinear combination of two modes with wavevectors fci and fc2 , and 
with Q{k\) = m/2 (blue), can, in the presence of a second component of the forcing at frequency n, 
interact with a mode with wavevector fcs (red), provided ks = fci + fc2 and U{k[i) = \m — n\. Waves 
driven by frequency n (green) do not enter the resonance condition. 



applies, for example, to two-frequency forcing with the third wave oscillating at the 
difference frequency: ^2(^3) = m ~ n and A = n. Other cases analogous to (2), such 
as n(/c3) = m + n, are discussed in [61], where the special significance of the differ- 
ence frequency case is explained. Note that in both cases (1) and (2), the temporal 
frequency ^^(fca) determines the angle 9 between the wave- vectors fci and fc2 via the 
dispersion relation, and therefore provides a possible selection mechanism for certain 
preferred angles appearing in the power spectrum associated with the wavepatterns. 
Selecting an angle of 0° , with 1 : 2 resonance in space and time (figure 12.1b ) , is a 
special case. 

The nonlinear interactions of the modes can be understood by considering reso- 
nant triad equations describing small amplitude standing wave patterns, which take 
the form 

Zl = Xzi + qiZ2W + (a|/;ip + b\z2\'^)zi H 

Z2 = Xz2 + qiziw + {a\z2\'^ + b\zi\^)z2 + ■■■ (2.2) 
w ^ lyw + q-iZiZ2 + • • • , 

where all coefficients are real, and the dot refers to timescales long compared to the 
forcing period. Here the quadratic coupling coefficients qj are 0(1) in the forcing 
in the 1 : 2 resonance case, and 0(|/„|) in the difference frequency case [49]. For 
other angles 9 between the wavevectors fci and fc2 we expect gi « 53 « because the 
temporal resonance condition for the triad of waves is not met. 

Since the m frequency is dominant, the zi and Z2 modes will be excited (A in- 
creases through zero) while the w mode is damped {v < Q in (|2.2p ). In this case, 
w can be eliminated via center manifold reduction [26] near the bifurcation point 
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(w ~ ^^^Tp'); resulting in the bifurcation problem 

ii = Azi + (a|zi|^ + 5|z2p)zi 

Z2^ Xz2 + ia\z2f +b\zif)z2 . (2.3) 

These equations describe the competition between standing waves separated by an 
angle 9, where b = b+^j^ explicitly includes the contribution from the slaved mode w, 
and hence depends on the angle between the two wavevectors fci and fc2- 

The contribution of the damped w mode is significant whenever gigs is non- 
negligible and the damping ji^j is not too great, and can be made more important 
by increasing qiq^ and/or by decreasing the damping ji^l. For instance, in the 1 : 2 
resonance case, for which qiq^ is 0(1), the damping can be decreased by taking n ~ 2m 
in (|2.ip . since the w mode, with frequency m, will then be driven subharmonically by 
the n component of the forcing (as well as harmonically by the m component). 

In the difference frequency case, the quadratic interactions rely on the presence of 
the n component of the forcing to allow the temporal resonance condition to be met. so 
qiq^ is 0(1 /„p). Thus the contribution of the damped w mode to b can be made more 
important in two ways: first, by increasing /„, or second, by parametrically driving 
the damped mode so that \v\ is decreased. This requires a third driving frequency 
p ~ 2\n — m| to be added to the forcing (|2.ip . In both of these instances, the relative 
phases of the components of the forcing matter, since the modes are being generated 
nonlinearly with a preferred phase. These ideas arc developed systematically in [49, 
61], where it is also shown that if there is an underlying Hamiltonian structure, then 
qiqs < for the 1 : 2 resonance, and qiqs > in the difference frequency case 
provided ji > m. Note that when gigs > (gigs < 0) then the fcs-mode makes a 
positive (negative) contribution to the growth rate of the mode Z2 when zi is present, 
and vice versa. 

In order to determine more precisely whether the resonant contribution to b is 
significant enough to lead to a qualitative change in the resulting pattern, it is useful 
to rescale the amplitudes Zi and Z2 in (|2.3p by a factor of l/-\/[a[. Then we obtain, 
for a < 0, the rescalcd equations 

ii = Azi-(|zi|2 + Be|z2nzi 

Z2 = Xz2~i\z2\^+Bg\zi\'')z2 , (2.4) 

where Bg = ^. Here the 6 subscript indicates that the cross-coupling coefficient 
between the fci and k2 modes depends on the angle 9 between them. 

The Be function has important consequences for the stability of regular patterns. 
As a simple example, note that stripes (|zi| = y/X, Z2 = 0) are stable if Bg > 1, while 
rhombs associated with a given angle 6 (|zi| = \z2\ = •\/A/(l -t- Bg)) are preferred 
if \Bg\ < 1. Moreover, if \Bg\ < 1 for any angle 9, then stripes will necessarily 
be unstable near onset. Since, by judicious choice of forcing frequencies we have at 
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least some ability to control both the magnitude and sign of Bg over some range of 
angles 0, we have a mechanism for enhancing or suppressing certain combinations of 
wavevectors in the resulting weakly nonlinear patterns. Alternatively, as suggested 
by [67] , if we choose forcing frequencies that lead to a large 1 : 2 resonant contribution 
at = (figure 12.1b ) , then the coefficient a can become large compared to &, which 
in turn will cause the rescaled cross-coupling coefficient Bg to be small over a broad 
range of angles away from = 0. (As ^ 0, it can be shown that Bg 2.) These 
two cases are contrasted in sections [7] and |8l with preliminary work described in [52] . 

Before continuing, we reiterate the constraint on using this analysis to design forc- 
ing functions that will stabilise a given pattern. In eliminating the damped mode w, 
we performed a center manifold reduction, which is valid provided that all non-neutral 
modes are linearly damped with decay rates that are bounded away from zero. The do- 
main of validity of the reduced equations depends on the extent to which w is damped 
near the bifurcation point. If that damping is very weak, then the reduced equations 
will only be quantitatively predictive in a correspondingly small neighborhood of the 
bifurcation point. We point this out since, from the discussion above, it is clear that 
the two ways of influencing the magnitude of fe, namely increasing the driving force fn 
to control (71(73, or driving the difference frequency to reduce the damping \v\, both 
can lead to situations where the center manifold reduction is no longer valid: the m 
and n modes could set in together, or the m and difference frequency modes could set 
in together. In either codimension-two analysis could be performed, but this 

is beyond the scope of this paper. In practice, the severity of this constraint can only 
be seen by comparing predictions for the amplitudes and stability of patterns with 
solutions of the problem at hand (which we do systematically in section [S]). 

3. The model PDE. In order to explore these issues in detail, we have devised 
a phenomenological PDE for which the leading nonlinear coefficients in bifurcation 
problems such as p.3p . can be calculated relatively easily, and which is reasonably 
simple to integrate numerically. The equation is: 

dU 

— = (/i + iuj)U +{a + ip)W^U + (7 + iS)\/^U 

+ QiU^ + Q2|C/|2 + C\U\^U + iRc{U)f{t), (3.1) 

where f{t) is a real 27r-periodic function, U{x,y,t) is a complex- valued function, with 
{x,y) € M^, and // < 0, w, a, /?, 7 and S are real parameters, and Qi = Qu- + iQu, 
Q2 = Qir + ^Qii and C = Cr + iCi are complex parameters. 

The way the forcing function enters the PDE was chosen so that the linearised 
problem reduces to the damped Mathieu equation (in much the same way that hydro- 
dynamic models of the Faraday instability reduce to this equation [3]). The PDE has 
the advantage that the dispersion relation can be controlled easily, and weakly nonlin- 
ear theory is relatively straightforward to compute. The linear terms are diagonal in 
Fourier space, so the PDE is readily amenable to the Exponential Time Differencing 
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numerical methods of [11]. In addition, the nonhncar terms are simple (they do not 
involve any derivatives), and so numerical solutions are relatively inexpensive. 

One of the special features of parametric systems is that, even though pattern 
selection is a nonlinear process, the position of the linear stability curves determines 
which resonant triad interactions are dominant. In turn, it is the resonant triad inter- 
actions rather than the details of the particular form of nonlinearity in the equation 
that drives the pattern selection process. For these reasons, the model PDE is a useful 
testing ground for results derived from symmetric bifurcation theory. 

The model PDE is similar to the complex Ginzburg-Landau equation - but we 
point out that U{x,y,t) is itself the pattern-forming field, and is not the amplitude 
of some other underlying pattern. With ^ < 0, all waves are damped in the absence 
of driving. Note also that the parametric forcing f{t) is explicitly a function of time 
so that we are resolving dynamics on the fast time-scale set by the periodic forcing. 
In contrast; other authors [8-10, 28] have investigated Ginzburg-Landau equations 
that describe the slow, large spatial-scale evolution of the amplitude of an otherwise 
spatially homogeneous oscillatory mode arising through Hopf bifurcation. Instead 
of resolving the fast oscillations of the subharmonic response to the time dependent 
forcing, a constant-coefficient U term is introduced into the equation, proportional 
to the amplitude of the parametric forcing [10]. With multi-frequency forcing, other 
terms, such as C7^, are also introduced [8,9,28]. 

There are important qualitative differences, of course, between the model PDE 
and the Faraday wave experiment. One difference is the role of the fc = mode. In 
the PDE, this mode is damped (fi < 0) and has a non-zero frequency uj; moreover 
the fc = mode can be nonlinearly excited. In the Faraday wave problem, owing 
to mass conservation, the fc = mode is neutral and cannot be excited, and this 
may have important consequences in the dynamics [12,39]. (In the Zhang- Vihals 
model [67] this requirement is met since all nonlinear terms appear with an overall 
spatial derivative that prevents the excitation of the fc = mode.) Another important 
difference between the model PDE and the Faraday wave experiment is that the 
dispersion relations have a different structure: in the model PDE, the frequency is 
a polynomial function of the square of the wavenumber, but the dispersion relation 
for Faraday waves is more complicated [3,4,35]. Nonetheless, the marginal stability 
curves of the model PDE, especially with multi-frequency forcing, are similar to those 
that are observed in large aspect ratio Faraday wave experiments. 

The theory developed by Porter, Topaz and Silber [49,61] applies in the weakly 
damped, weakly forced regime, and certain of their results also require that the un- 
damped problem have a Hamiltonian structure. This limit can be realized for our 
model p.ip by setting fi = a = j = Cr^O and Q2 = — 2Qi- In this case, the 
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Hamiltonian is: 



H{U,U)^ // [c^|C/p-/3|V[/p+5|v2c/|2 + /(t)(Re(t/))2 



R 



and U evolves according to 



iQiU^U + iQiU^U +^\Uf] dxdy, (3.2) 



dt SU ^ ' 



The region R corresponds to the domain of integration of the PDE (|3.ip . where we 
have assumed periodic boundary conditions apply on dR. In the examples presented 
below, some of the parameter choices are nearly Hamiltonian (sections [7] and [8]) and 
some are not (section ^ . 

4. Linear theory. The linear problem associated with (|3.ip takes the form of a 
damped Mathieu equation for each Fourier mode e**"'^. Specifically, if we set U {x, t) = 
e*'^'^(|)(t) + iq{t)) in (j3.ip linearized about [/ = 0, then we obtain 

p + 7?j+(o2+O/(t))p = 0, (4.1) 

where 

" 2 

7 = 2 (-^ + ak"^ - jk'^) , n = uj - pk"^ + Sk'^, il"^ = ^ + n'^ . (4.2) 

We use the method of [4] to solve this linear problem for multi-frequency forcing f{t), 
which determines the critical forcing amplitude. Further details can be found in 
Appendix [Xl 

We require that the damping 7 be positive for all fc. so /i < 0, 7 < and 
a > — ^47/i. However, we do not necessarily insist that it be monotonic, which would 
require a > 0. Non-monotonic damping is possible in the Faraday wave experiment 
in shallow layers: if the viscous boundary layers extend from top to bottom of the 
experimental container, modes with long wavelengths can be more heavily damped 
than short-wave modes, and indeed the first harmonic mode can be unstable at lower 
forcing that the subharmonic mode [34,44]. We use non-monotonic damping for the 
example in section [71 Likewise, we restrict ourselves to the parameter regime where 
the dispersion relation n(k) is non-negative and an increasing function of fc^; thus we 
assume oj > 0, 5 > and (3 < 0. 

We will use both 1 : 2 and difference frequency resonances to control how modes 
interact, so it is important to understand that these resonances impose further con- 
straints on the parameters in the dispersion relation. In the case of single frequency 
forcing /(<) = acos{t + (/)), we expect a subharmonic instability to set in first with 
increasing a at a critical wavenumber ki , which can be estimated by solving the equa- 
tion Cl{ki) = i. This estimate of ki is good provided the damping is not too large. 
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and that the damping is a (nearly) monotonic function of k. The wavenumber asso- 
ciated with the first harmonic instability is determined by solving Cl{k2) = 1. These 
calculations, together with simple trigonometry, determine that the angle 6 associated 
with a 1 : 2 resonant triad satisfies the equation ^2 = 2ki cos (|) , which has a solution 
provided k2 < 2fci. Choosing a scaling of x so that fci = 1, then we find that uj, /3 
and S are related to 6* by 

h{k = l)=uj- I3 + S=^ 



n l^k = 2cos l^-jj =cj-4/3cos^ l^-j + 16(5 cos^ = 1 . (4.3) 

In particular, the 1 : 2 resonance will be at 6* = 0° (fc2 = 2) if we choose 

Lo^^+AS, /3 =-1 + 5,5, (4.4) 

where S G [0, ^] ensures P < and hence a monotonic dispersion relation. The 1 : 2 
resonance moves to 6* = 90° {k2 = V2) if we choose 

uj^26, /3 = -i+3(5, (4.5) 

where we require 6 € [0, 

Next we consider the case of two-frequency forcing 

f{t) = F {am cos(mt + (/>„,) + a„ cos(rit -f- (j>n)) , (4.6) 

where m and n are coprime integers, (am,a„) are relative amplitudes scaled by an 
overall amplitude F, and {(j)m, <t>n) are phases. We focus on a resonant triad involving 
two critical modes with dominant frequency and a damped difference frequency 
mode with dominant frequency n — m. We will typically take m even and n odd with 
n > m and n — m < ^ (*-e-, ^ G (li§))- These conditions imply that the initial 
instability is expected to be harmonic, that the difference frequency mode decreases 
Bg at the resonance angle in the Hamiltonian limit, and that the difference frequency 
mode has a wavenumber (fcdiff) that is smaller than the critical wavenumber and hence 
it is not too strongly damped (at least in the case of monotonic dissipation) [49] . We 
estimate the critical wavenumber of instability as fc w fci , where we assume a scaling 
such that fci = 1. We then have that fci and fcdiff satisfy 

n{k = i) = uj-p + 6 = — 

^{k = fc"diff ) =^^~ Pkj-ff + Sk^.ff ^n-m , (4.7) 
which we can solve for lu and /?: 



^= 'V J ^+ ^fcdiff and P= \_\^ L+5{l + kl^). (4.; 

'^diff ^ '^AiS 
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Fig. 5.1. (a) Rhombic lattice with an angle 6 between the primary wavevectors. (b) Hexagonal 
superlattice, with an angle of 21.8° between the most closely spaced wavevectors. (c,d) 12-fold and 
14-fold quasilattices, up to 11th order and 7th order respectively [51]. See section^for a discussion 
of how these quasilattices are generated. 



Setting fcdiff = 2 cos (^), we can relate an angle in the power spectrum of the pattern 
with the wavenumber fcdiff of the damped mode associated with the resonant triad. 
For instance, for 6 = 150° (the complementary angle to 30°, which appears in the 
12-fold quasipatterns), we have fcdiff = V 2 — a/3. Alternatively, if 6* = 158.2° (the 
complementary angle to 21.8°, which appears in the simplest hexagonal superlattice 
patterns), we have fcdiff = 

5. Weakly nonlinear theory. Our weakly nonlinear calculations are aimed 
at determining the coefhcients of the leading nonlinear terms in finite-dimensional 
bifurcation problems associated with certain families of patterns in the plane, and the 
corresponding lattices of wavevectors. These finite-dimensional bifurcation problems 
allow us to rigorously compute the relative stability of various simple planforms, 
e.g., stripes v.s. squares, rhombs, hexagons, and also to calculate relative stability of 
superlattice patterns and hexagons, stripes, and certain rhomb patterns. Moreover 
these calculations lead to quantitative predictions of the amplitude of the standing 
wave patterns as a function of the distance A from the bifurcation point, where A = 
{F — Fc)/Fc and Fc is the critical value of the overall forcing amplitude. 

A simple example of such a reduction to a finite-dimensional problem was pre- 
sented in Section [2l where we described the bifurcation problem associated with a 
pair of standing waves oriented at an angle 6 relative to each other, where 9 £ (0, 7r/2] 
was bounded away from 7r/3. An example rhombic lattice is shown in figure ISTlT a'). 
In that case, the bifurcation problem consisted of a pair of amplitude equations given 
by (|2.4p . and, after appropriate scaling, there was just a single nonlinear coefficient 
Bg . The details of the (numerical) computation of this coefficient from the governing 
PDE p.ip is given in Appendix [A] In subsequent sections of this paper we present 
plots of Bg as a function of d for certain parameter sets and forcing functions f{t) 
used in our numerical simulations of the model PDE. 

As an additional, specific example we consider patterns associated with the hexag- 
onal superlattice that has been observed in several Faraday wave experiments [2,20, 
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33]. Equivariant bifurcation theory [23] was used to derive the form of bifurcation 
problem [17]. This bifurcation problem describes the long-time evolution of the twelve 
modes on the critical circle that are associated with patterns that tile a plane in 
hexagonal fashion. The critical Fourier modes associated with the 21.8° superlattice 
pattern are indicated in Figure [SHT b). We label these modes as follows: (21,23,2:5) 
are complex amplitudes associated with wavevectors separated by 120° and, together 
with their complex conjugates, they comprise the modes associated with one hexagon, 
while (z2, Z4, zg) and their complex conjugates are associated with a second hexagon, 
rotated by approximately 21.8° relative to the first. 

The form of the bifurcation problem associated with this hexagonal superlattice, 
to cubic order in the amplitudes, is: 

^ = Azi + Qzg^s - d^il' + Beo ([zal" + l^sl') + + Ssslzel' + 5821^21') 

^ = AZ3 + - (|Z3|' + Beo {\Z5\'' + \zi\^) + S22IZ6I' + S38IZ2I' + ^821^41') Z3 

^ = \Z5 + QziZ3 - (I25I' + Seo d^ll" + |Z3|2) + B22IZ2I' + ^381^41' + ^821^61') ^5, 

with similar equations for Z2, Z4 and zg; related by symmetry. We have shortened the 
labels of the angles to 22 instead of 21.8, etc. This label indicates the angle between 
pairs of modes, e.^., there is an angle of approximately 82° between z\ and Z2, while 
z\ and Z4 are separated by approximately 22°. The nonlinear coefficients Q and Bgo 
are computed from the governing PDEs by considering the problem of bifurcation 
on a simple hexagonal lattice involving a subset of the modes, while Bq is computed 
for an arbitrary B ^ 60° on a rhombic lattice. We can then read off i322, -833, etc., 
as required. The details are in Appendix |^ Note that we have assumed that the 
bifurcation to a stripe planform is supercritical so that we can rescale the amplitudes 
to make the self-coupling coefficient a = — 1. 

The standard planforms, namely stripes, rhombs (associated with each of the 
angles 22°, 38°, 82°), hexagons, and superlattice patterns, correspond to equilibrium 
solutions of these equations. The calculation of their linear stability proceeds in 
a standard fashion and is summarized in [17]. In fact, due to the presence of the 
quadratic term, all planforms bifurcate unstably, but because Q is typically very small 
for multi- frequency forcing and sufficiently weak damping [47], the planforms can be 
stabilized at small amplitude by secondary bifurcations. The superlattice patterns, 
which satisfy z\ = Z2 = ■ ■ ■ = zq, come in two varieties that bifurcate together and 
their relative stability is unresolved at cubic order. Specifically, hexagonal superlattice 
patterns are associated with Zj being real, while triangular superlattice patterns are 
of the form Zj ~ Re^^l^ , where R is the real amplitude [54]. Which of these two 
superlattice patterns is favored in a given situation depends on higher order terms in 
the bifurcation problem. We do not calculate the coefficients of these terms and merely 
lump the two types of superlattice patterns together in our bifurcation diagrams. 
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6. Numerical experiments: selecting superlattice patterns. Motivated 
by experimental observations [2,20,33] of superlattice patterns in the Faraday wave 
experiment with 6 : 7 forcing, in this section we use 6 : 7 forcing, with some additional 
forcing at frequency 2 to drive the difference frequency mode, in order to stabilise a 
21.8° hexagonal superlattice pattern. We carry out the linear and weakly nonlinear 
calculations to find parameter values for which hexagonal superlattice patterns are 
predicted to be stable, and solve the model PDE numerically to confirm these predic- 
tions. In addition, we make quantitative comparisons between the weakly nonlinear 
predictions and the numerical solutions of the PDEs, comparing the predicted am- 
plitudes and ranges of stability of the patterns. The agreement is not quantitative 
except at very small amplitude, and we develop an explanation for this at the end of 
the section. 

The PDE was solved numerically using the fourth-order Runge-Kutta Exponen- 
tial Time Differencing numerical method (ETD4RK) of [11]. This pseudospectral 
method solves the linear part exactly, excluding the parametric forcing term, which is 
included with the nonlinear terms. This allows the use of a timestep based on accu- 
racy requirements, rather than numerical stability limits. Timestepping takes place 
in spectral space, and the timestep is chosen to be one-twentieth of the shortest of 
the periods of the forcing function f{t), in order that the effect of the time-dependent 
forcing is fully resolved by the fourth-order method. The nonlinear terms are eval- 
uated using FFTW [21]. The resolution was relatively low for the examples in this 
section (up to 96 x 56 Fourier modes), but we used up to 1536^ Fourier modes for the 
largest quasipattern examples discussed below. At each timestep, the upper half of 
the Fourier spectrum was removed in order to dealias the cubic terms. 

We choose parameters so that the mode driven subharmonically by the 6 frequency 
has wavenumber k = 1: 0(fc = 1) = 3. The wavenumber associated with 21.8° is 
which we wish to correspond to the difference frequency (7— 6), so we set ^("^) = 1- 
This mode responds subharmonically to driving at frequency 2. Furthermore, the 
wavenumber associated with 81.8° is part of this superlattice pattern, and we can 
influence this if we set Cl{-^) = 6, in 1 : 2 resonance with the primary response. 
This yields a; = |, /3 = — | and 6 = 0. We choose damping coefficients /i — —0.1, 
a = 0.01 and 7 = —0.15, and nonlinear coefficients Qi = 2 + i, Q2 = 1 + 2i and 
C = —1 + 3Qi. Note that these parameters are not close to the Hamiltonian limit: in 
fact, we have chosen the nonlinear coefficients so that the 1 : 2 interaction reduces the 
cross-coupling coefficient in the range of angles close to 90°. 

The linear theory for this problem is shown in figure 16.11 confirming that modes 
with k close to 1 and are neutral and weakly damped respectively, with 6:7:2 
forcing (figure 16. IF ). For this example, we have set the phases of the two main 
components of the forcing equal to zero, and (after some experimentation) set the 
phase of the component that drives the difference mode equal to 240°. 

The cross-coupling coefficient Bg is shown in figure [§21 with 6 forcing only, there 
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Fig. 6.1. An example of the linear theory for one, two and three-frequency forcing, with dis- 
persion relation coefficients lu = 2/3, /3 = —7/3 and 5 = 0, and damping coefficients fi = —0.1, 
a = 0.01 and 7 = —0.1. (a,b) 6 forcing, with as = 1 and Fc = 0.83973. (c,d) 6 : 7 forc- 
ing, with (06,07) = (1,1), ((^6,07) = (0,0) and Fc = 0.80839. (e,f) 6:7:2 forcing, with 
(06,07,02) = (1,1,0.45), {<j>e,<f>7,<f>2) = (0,0,240°), Fc = 0.80975 and kc = 0.9910. (a,c,e) show 
neutral stability curves, whose minima define Fc and the critical wavenumber (close to 1 in all 
cases). Curves corresponding to the response to frequency 6 are blue, to frequency 7 are green, and 
to twice the difference frequency are red, with the corresponding driving frequency indicated. The 
minimum of the red curve in (e) is close to k = (b,d,e) show the real parts of the Floquet 

multipliers at the critical forcing. 



is a dip in the curve around 90° owing to 1 : 2 resonance at 82°. (Had we chosen 
parameters sufficiently close to a Hamiltonian limit, this feature would have been a 
peak rather than a dip [49].) With 6 : 7 forcing, the dip at 22°, corresponding to 
the difference frequency, is visible, even though this frequency is not forced directly. 
Finally, with 6:7:2 forcing, wc can control the depth of the dip at 22°. An 
additional feature at 43° is visible, corresponding to frequency 4. in 1 : 2 resonance 
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Fig. 6.2. Bg, with linear parameters as in figure [KT\ and nonlinear coefficients Q\ = 2 + i, 
Q2 = 1 + 2i, C = —1 + 30i. The blue curve is with 6 forcing only: the dip at 90° is because 
0/ 1 : 2 resonance with frequency 6. The green curve has 6 : 7 forcing: note the dip starting at 
22°, corresponding to the difference frequency, even though this frequency is not forced directly. 
The red curve has 6:7:2 forcing: note the pronounced dip at 22°, and the smaller dip around 
43°, corresponding to frequency 4 (which is in 1 : 2 resonance with frequency 2). Bqq is calculated 
separately, as described in Appendix [21 The relevant coefficients are B22 = 0.46, -B38 = 1.23, 
Beo = 1.11 and Bg2 = 0.63. 



Fig. 6.3. With 6:7:2 frequency forcing at (a,b) I.OO4 and (c) 1.02 times critical, we find 
hexagons, a superlattice pattern, and 22° rectangles. The resolution was 96 X 56 Fourier modes, and 
the domain is rectangular, 2\fl X 2^7 /Z critical wavelengths, big enough to fit two copies of the 
superlattice pattern. The grey-scale represents the real part ofU{x,y,t) with t equal to an integer 
multiple of 2tt. 

with frequency 2. 

Solutions of the PDE with 6:7:2 forcing arc shown in figure 16.31 confirming 
that hexagons, superhexagons and 22° rectangles are all stable solutions for different 
parameter values. 

In figure 16. 4[ we make quantitative comparison between the amplitudes and sta- 
bility of these patterns as numerical solutions of the PDE, and the values predicted 
by weakly nonlinear theory, for forcing up to 1.02 times critical. In addition, we show 
the results of one-dimensional simulations, which recover the unstable stripe pattern. 
In all cases, numerical solutions of the PDE agree with the weakly nonlinear pre- 
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Fig. 6.4. Amplitudes of hexagons (blue), the superiattice pattern (red), 22° rectangles (magenta) 
and stripes (black) as a function of \ for 6:7:2 forcing. Using weakly nonlinear theory, we denote 
stable solutions with solid lines and unstable solutions with dashed and dotted lines. We also denote 
numerically computed stable solutions of the PDEs by + for hexagons, X for the superiattice pattern, 
□ for rectangles, and o for (unstable) stripes. The amplitudes of the PDE solutions are computed 
by matching the time evolution of the k = 1 modes to the linear response functions, (b) is a detail 
of (a), showing that the agreement between weakly nonlinear theory and the PDEs improves very 
close to onset. 



diction, provided we are close enough to onset (we have confirmed this by plotting 
the data on a logarithmic scale). Secondary bifurcations, which delimit parameter 
intervals where the patterns are stable, are also recovered, although the agreement is 
only qualitatively correct. 

It is notable that the agreement between the amplitudes predicted by weakly 
nonlinear theory and measured from PDE simulations is not particularly good for the 
multi-mode patterns, when compared to the much better agreement in the case of 
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stripes. The reason for this lack of quantitative agreement for the complex patterns 
can be understood by going to higher order in the center manifold reduction that 
was performed to go from (|2.2p to (|2.3|) . We recall the framework: there are two 
weakly excited modes with amplitudes zi and Z2, and a damped mode w, which 
evolve according to: 

ii = Xzi + qiZ2W + (a|zip + b\z2\^)zi, 

Z2 = Az2 + qiziw + (a|z2|^ + b\zx\^)z2, (6.1) 
w = vw + q3ZiZ2, 

where all coefficients are real, and v < 0. We have discarded those higher order terms 
that do not play a role in the centre manifold reduction in order to emphasise the 
effect of the higher-order nonlinear terms that appear as a result of the reduction. We 
express w on the centre manifold as a power series in zi and Z2 , and perform a centre 
manifold reduction to find 

?3 , qsiqiQs -_{a + b)iy 

,'3 

which results in 

il = Azi + (a|zip + b\z2\')z^ + ^^-^3^.1.3-,^ r^y^y (j^^p ^ |^^|2) |^^|2^^ ^ 



ZlZ2^ ^ ( Zl + 22 jZl22H , (6.2) 

|2 , l\ \2^ , 9193(9193 - (a + 6);^) . 2 ,1 |2M |2. 



\ , ^ I |2 I ;;i |2x , 9193(9193 - (a + &)!^) /, |2 , I |2M ,2 , 
Z2 = AZ2 + (a|z2| +b\zi\)z2^ 5 (|2l| + \Z2\ ) \zi\ Z2 + ■ 



(6.3) 



where b = b — qiq^ / as before. In order for the w mode to influence the coupling 
constant Bg ~ b/ a, and so produce interesting patterns, the quadratic coefficients qi 
and 53 must be non-zero (the modes must be in three-wave resonance) and v must 
be small (the mode must be weakly damped). However, if v is small, the in the 
denominator of the quintic terms imply that these high-order terms will be important 
exactly where the most interesting patterns will be found. Indeed, the graphs of 
amplitude against A in figure are well fit by a quintic polynomial. 

In contrast, the stripe pattern involves damped modes at /s = and k ~ 2, and it 
can be seen from figure lHTlT f) that these modes are well damped, so there is no reason 
for quantitative agreement not to extend to larger values of A, and indeed it does. 

One can estimate the range of validity of the cubic truncation of (|6.3p when A 
and are both small. In the case of rhombs, the linear and cubic terms balance when 
|zip — |z2p — 0{Xiy). The quintic terms are thus smaller than the linear and cubic 
terms when A <C i^, which is what one would expect: the center manifold reduction is 
valid when all modes that are eliminated are heavily damped compared to the modes 
that are retained. 

Therefore, this codimension-one approach to finding interesting patterns has the 
smallest range of validity exactly where the patterns are most likely to be interest- 
ing. A proper treatment would require consideration of the codimension-two problem 
(A, ly) = (0, 0), which is beyond the scope of this paper. 
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A further complication in the Faraday wave situation (and in our model PDE) is 
that in order for the interaction between the primary harmonic modes (driven by the 
m forcing) and the weakly damped difference frequency modes (n — m) to take place, 
the subharmonic mode (n) must be present in the forcing function: the quadratic 
coefficients qi and qs increase with the strength of the subharmonic forcing /„, so 
the interaction is strongest when /„ is largest. This implies that the subharmonic 
mode is itself only weakly damped, and will be excited if /„ is increased beyond its 
critical value. Experimental evidence suggests that the codimension-two point (or 
bicritical point), where the primary harmonic (m) and subharmonic (n) modes are 
both neutral, is an organising centre for the dynamics [2,19,33]. (There has been 
relatively little progress on the theoretical understanding of the bicritical point, apart 
from a study in the case of a single frequency [65] and in a few particular cases for 
two- frequency excitation [47,48].) Therefore, the problem should really be treated as 
a codimension-three interaction between primary harmonic modes and subharmonic 
modes, as well as weakly damped harmonic modes. 

Notwithstanding these complications, it is clear that the idea that pattern se- 
lection is being influenced by three-wave coupling to weakly damped modes is fun- 
damentally correct, and the codimension-one approach, while having limited quanti- 
tative agreement with PDE simulations, is clearly providing the correct qualitative 
interpretation of the observed patterns. 

7. Numerical experiments: 12-fold quasipatterns. In the previous section, 
we demonstrated how to stabilise simple patterns by driving the difference frequency. 
In this section, we show how this mechanism can be used to predict parameter values 
for stable approximate 12-fold quasipatterns, and we demonstrate how well a periodic 
pattern in a large domain can approximate a quasipattcrn. 

In order to use triad interactions to encourage modes at 30°, we choose m = 4, 
n = 5 forcing: 4 : 5 forcing has been used in several experiments to produce 12-fold 
quasipatterns [1,19,33]. We set V.{k = 1) = 2 so that the subharmonic response 
to frequency 4 comes at wavenumber 1, and we require that a wavenumber involved 
in 150° mode interactions (fc^ = 2 — a/3) correspond to the difference frequency: 
Cl{k) = 1. One solution is w = 0.633975, (3 = -1.366025 and (5 = 0. Twelve- 
fold quasipatterns also require modes at 90° to be favoured, and for this choice of 
parameters, 0(fc = ^/2) is 3.37. Although this is not particularly close to 4, we can 
use 1 : 2 resonance (driving at frequency 8) to control the 90° interaction. The linear 
theory for these cases is in figure [7t1 with /i = —0.2, a = —0.2 and 7 = —0.15. Note 
that the damping is non-monotonic, and has a minimum at fc = 1. 

The resulting Bg curve with 4:5:8 forcing and with near-Hamiltonian choice of 
nonlinear coefficients (figure [7T2|) shows pronounced dips at 30° and 90° as required. In 
figure [7T3I we show that within a 12-amplitude cubic truncation, 12-fold quasipatterns 
are stable between 0.9995 and 1.0085 times critical. We have found stable approximate 
quasipatterns (marked by crosses on the figure, and discussed in more detail below) in 
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Fig. 7.1. Linear theory for the quasipattem example. The dispersion relation coefficients 
are w = 0.633975, /3 = —1.366025 and <5 = 0, and the damping coefficients are fi = —0.2, a = 
-0.2 and 7 = -0.15. (a,b) 4 forcing, with ai = 0.57358 and Fc = 1.04730. (c,d) 4 : 8 forcing, 
with (a4,a8) = (0.57358,1.6), {<p4,<f>a) = (0,0) and Fc = 0.95214. (e,f) 4:5:8 forcing, with 
(04,05,08) = (0.57358,0.81915, 1.6), {4>4,4>5,4>s) = (0,0,0), Fc = 0.93159 and fcc = 0.9798. (a,c,e) 
show neutral stability curves, whose minima define Fc and the critical wavenumber (close to 1 in all 
cases). Harmonic curves are blue and subharmonic curves are green, with the corresponding driving 
frequency indicated. (b,d,e) show the real parts of the Floquet multipliers. 



the same range. Squares are also stable above 1.0013 times critical. The agreement is 
good, better than the examples in the previous section (figure |6^ . since the difference 
frequency mode is still fairly well damped (compare the damping for k close to 0.378 
in figure [6TTF and k = 0.518 in figure Elf). 

7.1. Choice of domain size for approximate quasipatterns. Before pre- 
senting numerical solutions in large domains, we discuss the issue of choice of domain 
for providing accurate approximations to quasipatterns. 
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Fig. 7.2. Cross-coupling coefficient Bg for the parameters from figure \7A\ and nonlinear co- 
efficients Qi = 1 + i, Q2 = —2 + 2i, C = —1 + lOi. The blue curve is with 4 forcing only: the 
peak near 80° is because 0/ 1 ; 2 resonance with frequency 8. The magenta curve has 4 : 8 forcing, 
using the frequency 8 component to bring down the curve close to 90° . Bringing in the 5 frequency 
(red curve) give a pronounced dip at 30°, corresponding to the difference frequency (5 — 4), even 
though this mode is not driven directly. The relevant coefficients are B30 = —0.01, Bgo = 2.24 and 
B90 = -0.51. 

Just as hexagonal patterns can be approximated in rectangular domains [38], 
there are many ways of choosing periodic domains to allow accurate approximations 
of quasipatterns. Here we discuss three plausible approaches to choosing domains for 
12-fold quasipatterns. We show why one approach, based on Pythagorean triplets, 
docs not work at all, while two other approaches both work well. 

Reducible symmetry group representations with square periodic domains have 
been put forward as candidates for producing approximate quasipatterns [14,58]. In 
particular, Pythagorean triplets have been identified as of particular interest in this 
case [14]. Consider a pair of integers (p, q) with p > q > 0. Then (p^ — g^, 2pq,p'^ + q'^) 
forms a Pythagorean triplet {i.e., {p^ — q^Y + {2pq)^ = (p^ + q^)^), and the vectors 
fci = (liO) and fc2 = (2pq,p^ — (l^)/{p^ + Q^) have the same length. For example, 
with p = 7, g = 4, wc have fei = (1, 0) and k2 = (||, ||). If | is a continued fraction 
approximation of -^/S (table 17. ip , then the angle between these two vectors, namely 
tan^^ [{p^ ~ 9^)/ (2m)) > tends to 30° +0{l/q^) as the approximation to VS improves. 
Thus it might be thought that square (p^ + q^) x (p^ + domains might readily 
allow approximations to 12-fold quasipatterns. 

Unfortunately, these Pythagorean domains do not provide good approximations 
to quasipatterns. The reason is that the essential 60° coupling is not quite correct: 
consider the wavevector fei = {p^ + q^,0) (in units of the basic lattice vector), and, 
at 60° on cither side of fci there are the wavevectors fcs = (p^ — q^,2pq) and fen = 
(p^ — q'^, —2pq). In order for quadratic interactions between these three modes to 
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Fig. 7.3. Bifurcation diagram with 3-frequency 4:5:8 forcing based on a 12-amplitude 
cubic truncation, with parameter values as in figure [7^1 12-fold quasipatterns (red) are predicted 
to be stable up to 1.0085 times critical; squares (magenta) are stable from 1.0013 times critical 
and 30° rectangles (green) from 1.0092 times critical. The crosses show the amplitudes of stable 
approximate quasipattern solutions of the PDE, calculated in a domain 2\/T3 X 2-^13/3 critical 
wavelengths. 
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occur, we need fci = fca + fcn. However, in this case, fea + fen ~ {2{p^ — q^), 0) = 
ki + (p^ — 3q^,0), which is close to, but never equal to, fci: it can be shown that 
— 3q^ = 1 or —2, which is small is compared to |fci| — +q^- This small difference 
means that the important nonlinear interactions between these three waves generate 
erroneous long- wave modulations in square {p^ + q^) x {p^ + q^) domains, and if modes 
with wavenumber close to zero are not heavily damped, these long-wave modulations 
can dominate the pattern. 

A much better way of generating good approximations to 12-fold quasipatterns 
is to choose 2q x 2q domains, with vectors fei = (1,0) and k-2 = {p,q)/2q, with | 
drawn from table 17. Il Again, the angle between these vectors goes as 30° 4- 0{l/q^), 
and, though the wavcnumbcrs arc not quite equal, we have \k2\ ^ 1 -|- 0{l/q'^). 
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66.2 
-132.3 

66.2 
-132.3 

66.2 



p 
1 

7 

4 

19 
11 
26 
15 

11 
41 
97 
56 



2q X 2g 

64 

484 

900 

6724 

12544 



fci 
fci 
fci 

fci 
fci 



(1,0), fc2 
(1,0), fc2 
(1,0), fc2 
(1,0), fc2 
(1,0), fc2 
(1,0), fc2 



JL 1 

2g' 2 

- 

8 ' 2/ 
19 1 
22 ' 2, 
26 1 
30 ' 2, 
71 1 
82 ' 2, 



fci = (1,0), fc2 



• 97 1\ 
, 112 ' 2) 



V3 1 
2 ' 2 



0.498 
-1.001 

0.500 
-1.000 

0.500 



-16.3 
33.2 

-16.5 
33.1 

-16.5 



Table 7.2 

Domains that provide good approximations to 12-fold quasipatterns . The first column gives the 
rational number ^ that is a continued fraction approximation to \/3, drawn from table \ 7. 1\ The 
second and third columns give the area of a computational domain and two of the wavevectors that 
will make up an approximate quasipattem. The fourth and fifth columns demonstrate that the errors 
in the length of k2 and in the angle Z12 between fci and fc2 scale as 1/A, or equivalently, as q~^. 
The first set of rows refer to rectangular domains of size 

which allow superlattice patterns that approximate quasipatterns. In these examples, all wavevectors 
are the same length. The second set of rows refer to square domains of size 2q X 2q. In these 
domains, approximate quasipatterns have wavevectors that have two slightly different lengths. 



This approach works because with this choice of wavevectors, we do have the correct 
three-wave couphng: fca = {q,p)/2q, fcn = {q, —p)/2q and so fca -|- fcn = fei. 

A third possibility is to approximate the quasipattern as a superlattice pattern, 
using the 12-dimensional irreducible representations of Dq ix [17]. If we choose 
a = q and f3 = p — q (in the notation of [17, table 2]), with ^ drawn from table FTTl 
then fci = {{p - q)y^,3q - p), fc2 = {{2q - p)Vi,p), ks = {{p - 2q)V3,p), fcn = 
{{qV3,3q — 2p) (all these should be divided by their length, \J p^ — "ipq + 3g^). We 
have fca -|- fcn = fei as required, all wavevectors are the same length (which is an 
advantage over the second alternative), and the angle between fci and fe2 goes as 
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30° + 0(1/(7^). One disadvantage of this approach is that numerical solutions of the 
PDE must be carried out in 2 -^/p^ — 3pq + 3g^ x 2-^/ (p^ — 3pq + 3g^)/3 rectangular 
domains in order to take advantage of spectral numerical methods. These domains 
are big enough to contain two repeats of the pattern (as in figure 16.3b ) and so only 
half the Fourier coefficients are used. 

The last two methods are compared in table [7T2l The error between the approxi- 
mation and the 12-fold quasipattern is inversely proportional to the area of the domain 
in both cases. The hexagonal superlattice approximations have all wavevectors of the 
correct length, unlike the square approximation. For this reason, the superlattice case 
is better for computing bifurcation diagrams, as all modes bifurcate at the same value 
of the forcing. 

However, for similar domain areas, the wavevectors in the square case are about 
four times closer to 30° apart than the superlattice case. Moreover, the square case is 
more amenable to efficient use of fast Fourier transforms, since the number of modes 
can be chosen to be a power of two in each direction, while the superlattice case leads 
to more awkward choices of numbers of modes. 

For these reasons, we choose the rectangular 2\/T3 x 2-\/13/3 domain for the 
bifurcation diagram in figure [7?3l since the quantitative comparison between numerical 
simulation and analysis is easier when all modes have the same wavenumber and hence 
bifurcate at the same value of the forcing (see below). We use the more convenient 
square (8 x 8, 30 x 30 and 112 x 112) domains for the remaining PDE simulations 
described below. It would be interesting to see how sensitive the quasipattern is 
to the exact choice of domain size. It would also be interesting to find domains 
that provide particularly accurate approximations to 8-, 10-, 14-fold and higher order 
quasipatterns, but these issues are beyond the scope of this paper. 

One consequence of the wavevectors being of unequal length in the square case is 
that the two wavenumbers concerned have slightly different critical forcings. In the 
30 X 30 and 112 x 112 cases, this difference is negligible, but in the 8x8 case, the 
critical forcings for the two wavenumber are appreciably different. We therefore make 
a small adjustment to the domain size to make the two critical forcings the same, 
raised by a factor of 1.00066 above Fc- 

7.2. Numerical examples of 12-fold quasipatterns. Numerical solutions of 
the PDE ()3.ip with 4:5:8 forcing at 1.003 times critical are shown in figure [731 in 
periodic domains 8 x 8, 30 x 30 and 112 x 112 wavelengths with periodic boundary 
conditions and using up to 1536^ Fourier modes. The solutions were followed for at 
least 10 000 forcing periods in the largest domain. Most initial conditions resulted 
in square patterns, but minor adjustments to the Fourier amplitudes at an early 
stage of the calculation resulted in stable approximate 12-fold quasipatterns. Note, 
however, that the PDE solutions in figure [7^ were not constrained to chose exactly 
the wavevectors given in table 17.21 

The accuracy of the approximation improves with increasing domain size. The 
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(a) (b) 




Fig. 7.4. Numerical solutions of the PDE with 4:5:8 forcing at 1.003 times critical, (a) SxS 
domain, (b) 30 X 30 domain, (c) 112 X 112 domain (only a 60 X 30 section of the domain is shown). 
Parameter values are as in figure \7.S\ 

modes involved in the 30 x 30 example are (30, 0) and (26, 15) and their reflections, 
which are 29.98° apart, and differ in length by 0.05%. The larger 112 x 112 domain 
allows an improved approximation to the quasipattern: the important wavevectors 
are (112,0) and (97,56), which are 29.9987° apart and differ in length by 0.004%. 
The amplitudes of these modes differ by 1%. We discuss other ways of evaluating the 
improved approximation to quasipcriodicity in the next section. 

7.3. Fourier spectra of quasipatterns. All the numerical PDE solutions pre- 
sented here have been carried out in periodic domains, so these solutions are only 
approximations to quasipatterns. In the 12-fold examples (figures 17.41 and 18.3b be- 
low), the most important twelve wavevectors in the pattern arc not exactly 30° apart 
and do not have exactly the same length (sec table \T7^ . On the other hand, it is 



Design of paramctrically forced patterns and quasipatterns 



25 



possible to generate true quasipatterns using twelve modes with fc = 1 evenly spaced 
around the unit circle, but the asymptotic series in the weakly nonlinear approxi- 
mation are known to diverge [51]. In this section, we compare the detailed Fourier 
spectra of the approximate quasipatterns to see how these differ from the spectra of 
true quasipatterns. 

We make the comparison by computing the locations of modes generated by 
nonlinear interactions up to a certain order, in the three cases of approximate quasi- 
patterns in figure TTM as well as in a true quasipattern. To do this, we first define the 
order of a mode. 

Quadratic interactions between the twelve modes with wavevectors fci, . . . , fci2 
generate new modes with wavevectors 2fci, fci -|- A;2, fci + ks and so on. Nonlinear 
interactions of N of the twelve modes generate modes with wavevectors km'. 

12 

km='^mjkj, (7.1) 

where the m^-'s are non-negative integers adding up to N. We define \m\ = X)j=i 
In the case of a periodic domain, the set of all possible km defines a lattice of wavevec- 
tors (figure [STTb .b). while in the quasipattern case, the set of all possible km defines 
a quasilattice: examples of 12- and 14-fold quasilattices with |m| < 11 (|m| < 7 in 
the 14-fold case) are shown in figure ISTlT c.d). 

If a wavevector k is in the lattice or quasilattice, then k = km for an infinite 
number of choices of vector m, since increasing mi and my (say) by the same amount 
does not change km but increases |m| by 2. We define the order of wavevector k to 
be the smallest value of |m| with km = k. 

In figures [7751 and FTBl we compare the locations of wavevectors in the 8x8 (cyan), 
30 X 30 (blue) and 112 x 112 (red) examples from the simulation results shown in 
figure [7^ The 8x8 spectrum is only in figure FTSf a) (up to 7th order), and the 
30 x 30 spectra are only up to 26th order. The amplitudes of each Fourier mode 
is given by the size of the symbol (on a logarithmic scale): the largest symbols are 
the modes with the largest amplitudes, and the modes with the smallest symbols 
have amplitudes 10"'^^ times smaller. Modes with amplitudes smaller than this are not 
plotted as these are in the realm of round-off error (see below) . 

The success with which the periodic patterns approximate a true quasipattern 
can be judged by the locations of the important modes in the pattern, as compared 
to the locations of modes on the quasilattice, up to a given order. Up to 7th order 
(figure [73k ). the 30 x 30 and 112 x 112 modes (blue and red) overlay each other almost 
exactly, and correspond well with the quasilattice modes (small black dots at the 
centre of the red markers). However, the 8x8 modes (cyan) deviate substantially from 
the correct locations, and we conclude that the periodic pattern in an 8 x 8 domain 
is a poor approximation to a quasipattern (in spite of appearances in figure 17^^ ). 

At 11th order (figure [775b ). the agreement between 30 x 30, 112 x 112 and the true 
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Fig. 7.5. Fourier spectra o/ 8 X 8 (cyan), 30 X 30 (blue), 112 X 112 (red) approximate 12-fold 
quasipatterns, as well as the 12-fold quasilattice (small black dots), up to (a) 7th, (b) 11th and 
(c) 1 5th order. The size of the coloured marker indicates the amplitude of the corresponding Fourier 
mode on a logarithmic scale, with the smallest markers having 10^^ times smaller amplitude than 
the largest. Only Fourier modes close to the unit circle (0.8 < |fc| < 1.2) are shown, with wavevector 
angles —10° < 6 < 40°. The horizontal line k = 1 represents the unit circle, while the curved line 
represents fcj; = 1. Wavevectors that come closest to the unit circle up to a particular order are 
labelled with X . 
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Fig. 7.6. Fourier spectra of 30x30 (blue) and 112x112 (red) approximate 12-fold quasipatterns, 
as well as the 12-fold quasilattice (small black dots), up to (a) 26th, (b) 39th, (c) 56th and (d) 94th 
order. Here we show wavevectors within 10° of {1,0), and with 0.85 < |fc| < 1.15. 



quasipattern is still good, while at 15th order (figure [73b ). the 30 x 30 modes deviate 
noticeably from the 112 x 112 and quasipattern modes. This is more pronounced at 
26th order (figure ITTBb ) . 

At 15th order (figure [TTSh ) . the agreement between 112 x 112 and the true quasi- 
pattern is excellent: every red marker has a black dot at its centre. The agreement 
is still very good at 26th order (figure 17.6b ) : the black dots are not quite in the cen- 
tres of the smallest red markers, and some black dots do not have corresponding red 
markers. This means that those modes, present in the true quasipattern at this order, 
have amplitudes in the 112 x 112 approximation that are too small to be plotted. The 
situation at 39th order (figurc lTIBb ) is similar, but it isn't until 94th order (figurc lTTBt l) 
that the true quasipattern modes miss the centres of the red markers entirely. 
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Fig. 7.7. Fourier spectra of the 112 X 112 (red) approximate 12-fold quasipattern, as well as 
the 12-fold quasilattice (small black dots). All modes with amplitude greater than 10""'^ times the 
maximum amplitude are plotted; the quasilattice is plotted up to 56th order. The competing effects 
of the 12- fold order and the square lattice in are apparent. 

In fact, a curious situation arises at 56th and 94th order: it is apparent that modes 
close to fc = (1, 0), and in particular modes on the line kx = 1, have amplitudes that 
are higher than might be expected, since usually the amplitudes of modes decreases 
with their order. These modes are discussed in more detail below. However, it should 
be noted that plotting only modes up to a certain order masks the effect of the 
underlying lattice in the numerical solutions. In figure 17. 7|, we show all modes in 
the 112 X 112 down to round-off error, and quasilattice modes up to 56th order. 
The underlying square numerical lattice is clearly seen in the strings of red markers 
emanating from each large-amplitude mode. These give an impression that the large- 
amplitude modes could be considered to be clusters of modes in Fourier space. 

In figure FTSl we show the range of amplitudes of the modes as a function of order, 
for the 8 X 8, 30 X 30 and 112 x 112 approximate quasipatterns. We note that the 
8x8 example deviates significantly from the other two at all orders, while the 30 x 30 
and 112 x 112 examples arc in fairly good agreement until 19th order or so, another 
indication of how the approximation improves with the domain size. We note that 
direct comparisons between the 8x8 case and the other two cases are not strictly 
valid, as the domain size in the 8x8 case had to be adjusted slightly to allow for the 
different wavenumbers in the pattern, as discussed above. 

Figure rTST b) shows the spectrum of the 112 x 112 example at all orders up to 
100. The range of amplitudes at a given order decays exponentially with order up 
to about 23rd order, but then levels off at the level of the round-off error: about 
10~^^. On top of this pattern, there is a peak in amplitude at 15th order, and broad 
peaks around 41st, 56th and 97th order. We attribute these to the presence of weakly 
damped modes, with k close to unity, that appear around these orders. With weak 
damping, these modes will amplify the numerical noise in the PDE solutions, and 
nonlinear coupling implies that the amplified noise will feed in to nearby modes. 
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Fig. 7.8. Fourier spectra o/8x8 (cyan), 30x30 (blue) and 112 X 112 (red) approximate 12-fold 
quasipatterns, as a function of order, (a) up to 25th order, (b) up to 100th order (112 X 112 only). 
The vertical lines indicate the range of amplitudes of modes at that order, scaled to the maximum 
amplitude at order 1. 



This is shown in figure [7T91 where we plot \\km\ — 1| as a function of order |m| 
for the different examples. We note the marked drop in this quantity in particular 
at 15th and 39th order. While the drops in ||fc„i| — l| do not line up exactly with 
the peaks in amplitude, we suspect that it is this marked change in the damping of 
modes appearing at these orders that is responsible for the amplification of noise at 
around the same order. 

Of course, for the numerical patterns in periodic domains, there is a lower bound 
to 1 1 fern I — 1 1 that does not depend on |m| (this lower bound occurs within figure [7^ , 
while for the true 12-fold quasipattern, there is no lower bound [51]. If \km\ ^ 1, then 



Ifc. 



1 > 



K 



m ■ 



(7.2) 



where K is an order-one constant. For some particular choices of wavevectors, this 
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Fig. 7.9. Minimum value of \\km\ — l[, as a function of order \m\, for 8x8 (cyan), 30 X 
30 (blue), 112 X 112 (red) approximate 12-fold quasipatterns, and for a true 12-fold quasipattern 
(black) [51]. There are drops in this minimum quantity at 4th, 7th, 11th, 15th, 26th, 39th and 56th 
order in all cases. Modes that are responsible for these drops are identified in figures \7. 5\ and \7. 6\ 
In the 30 X 30 case, the Pythagorean mode (0.6,0.8) (which is of unit length) is generated at 26th 
order. 



asymptotic limit is achieved: for example, if 

km = pk4 + iq- l)kg + {q + l)fcii = - VSq), (7.3) 

wherep and q are integers, then \m\ ~ p+^q and |fcinp = l + \/3<z)^- When | is a 
continued fraction approximation to V3 (table FTT]) . then ||fcm,| — l| ~ ^{p— V^q)^ ^ 
-p-, so \\km\ — l| < , where K2 and K' are order-one constants [51]. These 
particular vectors are not always the closest ones that can be found at a given order, 
but they demonstrate that modes approach the unit circle arbitrarily closely (and so 
are arbitrarily weakly damped) as the order of the mode increases. For the fractions 
listed in table 17. 11 we have |m| = p-\-2q = 4, 11, 15, 41, 56, 153, . . . , but in fact there 
are step decreases in the minimum of \\km\ — l| at |m| = 4, 7, 11, 15, 26, 39, 56, 94. 
This is the reason for the choices of orders in figures 17.51 and 17.61 

In summary, having identified which sizes result in the most accurate approxima- 
tions to 12-fold quasipatterns, we find that in the largest example (112 x 112), the 
locations of the Fourier modes of the approximation deviate significantly from those 
of the true quasipattern only beyond 26th order. At this point, the amplitudes of 
the modes have reached the level of numerical round-off, so going any larger than 
112 X 112 would not lead to any significant improvement in the approximation to a 
true quasipattern for these parameter values. The small divisors reveal themselves by 
amplifying the amplitudes of the Fourier modes at (or close to) the order at which the 
small divisor appears, but they do not appear to cause the amplitudes of the Fourier 
modes to become excessively large, at least at this level of forcing, and up to the 
maximum order (153) available in this domain. 
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Fig. 8.1. Cross- coupling coefficient Bg for single frequency forcing, with 1 : 2 resonance in 
space and time: lu = ^, 13 = S = 0, fj. = -0.005, a = 0.001, 7 = 0, Qi = 3 + 4i, Q2 = -6 + 8i, 
C = —1 + lOi, Fc = 0.024002 and kc = 0.9999. Note that Bg drops away sharply as increases, 
and is close to zero for 9 > 30°. The relevant coefficients are B30 = 0.088, Bqq = 0.014 and 
Bgo = 0.010. 



8. Turbulent crystals: quasipatterns using 1 : 2 resonance. In order to 
have 1 : 2 resonance in space and time with single frequency forcing (m = 1), we 
impose = ^ and 51(2) = 1, which leads to = -1+4(5 and [3 = — i + 55. 

We choose 5 — 0, small values for the damping coefficients fi, a and 7, and order one 
values for the nonlinear coefficients close to the Hamiltonian limit. The resulting cross- 
coupling curve is shown in figure 18.11 as explained in section [2l the 1 : 2 resonance 
in space and time has enhanced the self-coupling coefficient (by about four orders of 
magnitude compared to the previous cases), and so the cross-coupling coefficient Bg 
drops away sharply, and is close to zero for 6 > 30°. 

Within the restrictions of a 12-mode expansion, 12-fold quasipatterns are stable 
(figure[S21)- Indeed, at 1.1 times the critical forcing, in a 30 x 30 domain, the numerical 
solution of the PDE with random initial conditions is a stable 12-fold approximate 
quasipattern f figure [575b ). As above, the primary modes that make up the pattern 
are (30,0) and (26,15) and their reflections, in units of basic lattice vectors. The 
amplitudes of the 12 modes differ by 0.5%. The initial condition was not in any 
invariant subspace, and the PDE was integrated for 160 000 periods of the forcing. The 
agreement between the weakly nonlinear predictions and the computed amplitudes is 
not good (figure [H21) J which we expect since the k = 2 mode is weakly damped (see 
discussion above). 

However, there is no feature in the cross-coupling curve (figure 18. ip to indicate 
that modes at 30° should enjoy a special status. When the forcing is changed to 
1.3 times critical with this 12-fold quasipattern as the initial condition, we find that 
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Fig. 8.2. Bifurcation diagram showing the weakly nonlinear predicted amplitude (solid line) 
and amplitudes of approximate quasipattern solutions of the PDEs in an 8 X 8 domain (crosses). 

it is unstable, and is replaced (after a transient of 70 000 periods) by a stable approx- 
imate 14- fold quasipattern (figures [H3t> and 18.41 and animation online). In this case, 
the 14 modes are (30,0), (27,13), (19,23), (7,29) and theh reflections, differing in 
length by 0.5% and having angles within 1.5° of 360°/14. The amplitudes differ by 
about 10%. Both the 12-fold and 14-fold quasipatterns are stable at 1.2 times critical, 
in a 30 X 30 domain. 

More complex quasipatterns are also possible: calculations done in larger 60 x 60 
and 90 x 90 domains at 1.3 times critical, starting with random initial conditions, 
both yield an approximate 20-fold quasipattern (figure [HSt). However, this is not a 
particularly accurate approximation to a 20-fold quasipattern: in the 60 x 60 case, the 
20 modes are (60, 0), (57, 18), (48, 36), (37, 47), (21, 56), (1, 60), (-18, 57), (-33, 50), 
(-47, 37), (-56, 21) and their 180° rotations, which differ in length by 0.4% and which 
have angles within 2° of 360°/20. The amplitudes of the 20 modes differ by up to 40%. 
However, figure [8T5l shows an examination of the Fourier spectrum of figure ISTST c): 
there are 20 peaks close to /c = 1, of similar amplitudes and arranged roughly evenly 
around the unit circle. The 90 x 90 example is similar. We speculate that in other 
domains and at other forcings, 16-fold and 18-fold approximate quasipatterns could 
also be observed. 

The cross-coupling Bg curve suggests that modes that are more than about 30° 
apart do not influence each other (at least in an amplitude equation truncated at cubic 
order). This suggests that patterns containing many modes at essentially arbitrary an- 
gles might be expected - such patterns have been termed turbulent crystals by Newell 
and Pomeau [45]. It is hard to see why one quasipattern should be favoured over 
another, on indeed, why quasipatterns (with wavevectors evenly distributed around 
the fc = 1 circle) should be favoured over more complex patterns. 
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Fig. 8.3. For the parameter values from figure [gTH we find three different approximate quasi- 
patterns depending on the amplitude of the forcing and the size of the domain, (a) 1.1 times critical, 
30 X 30 domain: a 12-fold quasipattern. (b) 1.3 times critical, 30 X 30 domain: a 14-fold quasipat- 
tem. (c) 1.3 times critical, 60x 60 domain: a 20-fold quasipattern (only 2 of the domain is shown). 
An animation of the transition from (a) to (b), also showing details of the Fourier spectrum, can be 
found online. 

We have not attempted to compare the locations of the high-order modes in these 
14- and 20-fold approximate quasipatterns with those in the true quasipatterns, as 
these higher quasilattices are very densely populated ffigure l^TTH ). and the approxi- 
mate solutions are not close enough to the true quasipatterns. 
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Fig. 8.4. Amplitudes of Fourier modes as a function of time, at 1.3 times critical m a 30 X 30 
domain. The initial condition was the 12-fold quasipattern from 1.1 times critical (red). This is 
unstable and, after an extended transient of 70 000 periods, it is replaced by the 14-fold quasipattern 
(green). Amplitudes of other Fourier modes close to k = I are shown in black. An animation of 
this transition can be found online. 
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Fig. 8.5. Amplitudes of Fourier modes of the 20-fold quasipattern from figure IKW c). as a 
function of wavevector orientation, for 0.95 < |fc| < 1.05, showing twenty peaks roughly evenly 
distributed. 

As an aside, there is an interesting connection that can be made between the 
Fourier spectra of these high-order quasipatterns and the fractal dynamics of the 
complex ODE (PC^/df?' = ~C"^^; where n is an even integer [25]. Solutions of this 
equation lie on multiply branched Ricmann surfaces, with branch points occurring 
densely at points in a quasilatticc of order n. 

The rotational degeneracy of the plane (in the absence of boundary conditions) 
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implies that any mode with |fc| = 1 is linearly excited, and it is an open question as 
to why patterns and quasipatterns, with a finite number of modes evenly distributed 
around the unit circle, should be the preferred patterns close to onset in so many exam- 
ples of pattern-forming systems. (Of course, pattern-forming systems are by definition 
those that produce regular patterns close to onset!) Newell and Pomeau [45] wrote 
down the evolution equation of N modes, ignoring quadratic terms and truncating at 
cubic order: 

, . N 

^ = XA,-Y^BoJA,fA,, (8.1) 

k=l 

where Aj is the complex amplitude of mode j, A is the growth rate, and 9jk is the 
angle between the wavevectors of modes j and k. Truncated in this way, the phases of 
the amplitudes do not enter the dynamics: Newell and Pomeau [45] attribute this to 
the translation symmetry of the underlying problem, but Melbourne [40] points out 
that this phase-invariance is in fact a normal form symmetry and hence not exact. 
This system of ODEs is variational, and evolves to minimise a free energy 

N N 

^ = -A^l^,f + - ^ BeJA,\Mk\'. (8.2) 

Newell and Pomeau [45] claim that in the case that the state that minimizes the free 
energy has many modes, with the magnitudes of all amplitudes equal but with 
arbitrary phases, then the pattern will resemble a spatially random field (a turbulent 
crystal). This case will be realised when Bg < 1 over a wide range of 9 (as in figure [8Tl 
for example). We therefore propose the quasipattern solutions in figure [8731 (12-, 14- 
and 20-fold quasipatterns) as examples of turbulent crystals. 

It is interesting to consider what determines the number of modes in these turbu- 
lent crystals. One argument, explored in more detail by Zhang and Vinals [68], is that 
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the preferred pattern at onset should be the global minimum of the free energy T . 
Assuming equal amplitudes, the free energy of an A^-fold pattern or quasipattern de- 
pends on Be evaluated at 360°/iV, 720° /N etc. Figure [HSl shows J- derived from the 
Be curves in figures [STT] and [7T2l with single and multi- frequency forcing respectively. 

In the single-frequency case (blue curve in figure [8?6)) . where Bg is close to zero for 
9 > 30° ffigure lSTTj) . there is a broad minimum around 14-fold or 16-fold quasipatterns. 
However, 12-fold, 14- fold and 20- fold examples were found (figure [8?3|) . depending on 
the forcing strength and domain size. No doubt other patterns could also be found 
with more exploration. 

In the multi- frequency case (red curve in figure 18. 6p , where Bg has pronounced 
dips at 30° and 90° ffigure [772]) . there are local minima at iV = 4 (squares), = 8 
and N = 12. In numerical experiments, most initial conditions found squares, though 
12-fold quasipatterns were also stable (figures 17731 and rr4|) . We have not looked for 
8-fold quasipatterns. 

These results suggest that these free energy arguments provide a useful qualitative 
tool for understanding pattern selection, but reality is often more complicated than 
the arguments might suggest. 

9. Conclusions. We have introduced a new model PDE (j3.1|) for investigating 
pattern formation and pattern design in parametrically forced systems. The PDE 
is intended to play the same role for the Faraday wave experiment that the Swift- 
Hohenberg equation [60] plays for convection: while the model cannot be derived from 
the fluid mechanics, it has qualitatively correct linear behaviour and the right type 
of nonlinear interactions in order to provide useful illumination of the processes that 
are going on. The model produces superlattice patterns (section [5]) and quasipatterns 
(sections [7] and [8]) in response to single and multi-frequency forcing, for the same 
reasons that these complex patterns are found in the Faraday wave experiment - 
confirming that the mechanisms are generic. The ease of calculating weakly nonlinear 
coefficients and of computing large-scale numerical solutions has allowed a quantitative 
exploration of the agreement between the theoretical understanding of the pattern 
selection mechanism and the patterns that are actually found. 

Of course, the model does not capture every detail of the physics of the Faraday 
wave experiment. In particular, the dispersion relations have different structures, and 
the model does not include the mean fiow effects that are important in nearly inviscid 
Faraday experiments [31]. The latter could be addressed by coupling the model to a 
conserved quantity or to a mean-flow equation (as in [24,63]) or by taking the negative 
Laplacian of the right-hand side of the PDE (as in [12, 14]). 

The Zhang- Vinals [67] equations do not have these drawbacks: these are derived 
from the Navier-Stokes equations in the limit of infinite depth and zero viscosity. One 
might ask what is gained by looking at a simpler PDE that is even further from the 
physics. There are two advantages of the new model: one is that it is very simple: 
the dispersion relation can be controlled easily for studying any resonant interaction 
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or response to multi-frequeney forcing; in addition, the weakly nonlinear theory can 
be computed very easily. A second advantage is that is is very well suited to the 
use of efficient numerical methods such as Exponential Time Differencing [11]: the 
linear terms are diagonal in Fourier space, and the nonlinear terms do not involve any 
derivatives. In contrast, the Zhang- Vinals equations are considerably more compli- 
cated and the weakly nonlinear computations are more involved. Numerical solutions 
are also more time-consuming, as the linear term is not diagonal in spectral space, 
and most of the nonlinear terms involve products of derivatives, resulting in more 
Fourier transforms for their evaluation. As a result of the relatively low cost of the 
calculations, we have been able to follow branches of solutions in detail, and to go to 
much larger domains and for much longer times than previous calculations. Of course, 
in the end it would be desirable to work directly with the Navier-Stokes equations, 
but for these, the weakly nonlinear theory is very challenging [57] and there are as 
yet no large-scale numerical simulations. 

Like the Zhang- Vihals equations, the new model includes explicit time depen- 
dence. In contrast; other approaches, based on developing a description of the 
slow evolution of the amplitude of an underlying pattern, use the Ginzburg-Landau 
equation with additional complex conjugate terms to capture the effect of the time- 
dependent forcing [8-10,28]. As a result, any complex patterns that are found must 
be interpreted in terms of slow, long-wavelength amplitude modulations of an under- 
lying pattern, which complicates any effort to make quantitative comparison between 
theoretical ideas and the behaviour of the real system. 

The significance of three-wave coupling to weakly damped modes and its role in 
pattern selection has long been recognised [19,41,45,68]. We have put this idea to a 
quantitative test by using it to choose forcing functions that stabilise a desired pattern 
in large domain calculations. However, the codimension-one approach to this idea, 
where the weakly damped modes are slaved to the pattern-forming modes, does not 
provide quantitative predictions of amplitudes of patterns, and of parameter regimes 
where the desired patterns should be stable, except for very close to onset. The 
reason is that computation of the cross-coupling coefficient Bg is only valid when 
all modes are strongly damped compared to the pattern-forming modes. This poses 
difficulties because the most interesting patterns occur where the pattern-forming 
modes are coupled to weakly damped modes, and this is where the theory used to 
calculate properties of these patterns is of limited validity. As a result, parameters 
had to be chosen very close to onset in order to find stable numerical examples of 
the desired patterns in parameter regimes where they were predicted to be stable. A 
codimension-two approach would extend the range of validity of the theory, and will 
be the subject of future work. 

We have investigated two mechanisms for the formation of quasipatterns. One 
mechanism uses three- wave interactions involving a damped mode associated with the 
difference of the two frequencies in the forcing to select a particular angle (30° in the 
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example presented here). Using different primary frequencies, or altering the disper- 
sion relation, would allow other angles, or combinations of angles, to be selected. The 
advantage is that a forcing function can be designed to produce a particular pattern: 
the mechanism is quite selective, and requires some fine-tuning of the parameters. 

The second mechanism uses 1 : 2 resonance in space and time to magnify the 
self-interaction coefficient and thereby, on rescaling, diminish the cross-coupling co- 
efficient Bg for angles greater than about 30°. This can lead to the formation of 
turbulent crystals [45]. The mechanism is robust (the patterns are found well above 
onset), and requires only single frequency forcing. A dispersion relation that supports 
1 : 2 resonance in space and time is needed. Within this framework, an inherent com- 
plication is that it is not clear why regular 8, 10, 12 or 14-fold quasipatterns, or indeed 
any other combination of modes, should be favoured. Indeed we have found that 12- 
, 14- and 20-fold approximate quasipatterns can be stabilised by altering the level 
of the forcing or the domain size, without changing other parameters, and we have 
reported the transition between two different types of quasipattern. The Lyapunov 
function approach [67] cannot make this distinction, and would predict that 14- or 
16- fold quasipatterns should be found at onset for these parameter values (figure [8^ . 
It remains an open question as to why one turbulent crystal should be favoured over 
another. 

The existence of 14- fold (and higher) quasipatterns has been suggested before [51, 
59,61,67], but we have presented here the first examples of spontaneously formed 
14-fold and 20-fold approximate quasipatterns that are stable solution of a PDE (fig- 
ure [8]3]), with preliminary results in [52]. Examples where 14- fold symmetry is im- 
posed externally have been reported in optical experiments [46]. The Fourier spectra 
of 12-fold and 14- fold quasipatterns are both dense (figure lSTh .d). but those of 14- fold 
quasipatterns are much denser, owing to the difference between quadratic and cubic 
irrational numbers [51]. This difference may have profound consequences for their 
mathematical treatment. 

We have identified what domain sizes result in the most accurate approxima- 
tions to 12-fold quasipatterns, based on square and on hexagonal domains ftable [7?2|) . 
and produced exceptionally clean examples of approximate quasipatterns in relatively 
large computational domains. Comparing the Fourier spectra of the approximate 
quasipatterns as a function of domain size, we have identified at what order the lo- 
cations of Fourier modes in the approximate quasipatterns deviate from those of the 
true quasipatterns. In the largest example (112 x 112), the locations of the Fourier 
modes deviate significantly only beyond 26th order (figure 17.6b ) , at which point the 
amplitudes of the modes have reached the level of numerical round-off (figure 17. 8p . 
This suggests that going any larger than 112 x 112 would not lead to any significant 
improvement in the approximation to a true quasipattern, at least for these parameter 
values. 

We have compared the amplitudes of the Fourier modes of the approximate quasi- 
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patterns and the leading order weakly nonlinear prediction f figures [7.31 and l8.2p . and 
found quantitative agreement very close to onset, but only qualitative agreement at 
larger amplitude, which is what would be expected from the problem of eliminating 
weakly damped modes, as discussed above. We have not extended this comparison to 
higher order since the weakly nonlinear calculations are too difficult for this paramet- 
rically forced problem. An extension of this work would be to devise a PDE without 
parametric forcing that also produces stable quasipattern solutions: this would allow 
high order weakly nonlinear calculations (as in [51]) and very large domain numeri- 
cal solutions, and so allow a comparison between computed mode amplitudes (as a 
function of order) and the weakly nonlinear theory. Standard weakly nonlinear the- 
ory produces amplitudes that diverge at high order because of the presence of small 
divisors [51], while the PDE solutions have amplitudes that decay exponentially with 
order - although the small divisors in this problem do make themselves felt by am- 
plifying the magnitudes of the Fourier modes at (or close to) the order at which the 
small divisor appears (figure 17. 8p . Such a PDE could be based on (for example) the 
Swift-Hohenberg equation [22,37,42], but the Swift-Hohenberg equation itself does 
not allow the weakly damped modes that are necessary to stabilise quasipatterns. 

Other numerical studies of quasipatterns as solutions of a PDE have not made 
a systematic study of the effect of domain size. Zhang and Vinals [67, 69] report 
approximate 8-fold quasipatterns in a 64 x 64 domain in their quasipotential model 
of the nearly inviscid Navier-Stokes equations, for parameter values close to the 1 : 2 
resonance in space and time. The modes involved were separated by 41°, 42°, 48° 
and 49° [69], so the approximation was not particularly accurate; our careful choice of 
domain size allowed much closer approximation. Miiller [42] developed a model based 
on two coupled Swift-Hohenberg equations, with parameters chosen so that the two 
unstable modes had wavenumbers that would favour 8-fold or 12-fold quasipatterns. 
Numerical simulations in a 10 x 10 domain in the second case found approximate 
12-fold quasipatterns. The modes involved are not stated, but we estimate them to 
be fci — (10,-1), k2 ~ (9,4), fca = (6,8) in units of the fundamental lattice vector. 
These have lengths 10.05, 9.85 and 10.00 respectively, and they are separated by 
28.48° and 29.17°. With fen = (4, -9), we have fes + fen = fei, so the 60° resonance 
condition is satisfied. Frisch and Sonnino [22] present a similar model and report 10- 
fold quasipatterns. Lifshitz and Petrich [37] found a 12-fold approximate quasipattern 
in a roughly 30 x 30 domain, in a model based on a single Swift-Hohenberg equation 
with a degenerate double minimum in its marginal stability curve. The modes involved 
appear to be the same as those in the 30 x 30 examples discussed in section [T] 

While we have not discussed the possibility of long-wave instabilities of quasipat- 
terns, the Fourier spectra of the 112 x 112 example (figures FTBl d) and l7.7p suggests 
that long-wave modes that are close to the primary wavevectors in a direction tangent 
to the critical circle are forced by high-order nonlinear interactions. This is also ap- 
parent from the locations of modes responsible for the small divisors [51]. Therefore, 
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any treatment of the long- wave stability of quasipatterns should take into account the 
presence of these modes. This is a delicate question. The only study of the sideband 
instabilities of quasipatterns [18] focusses on instabilities associated with modes that 
are perpendicular to the unit circle, using coupled Ginzburg-Landau equations for 
each of the primary mode directions in the quasipattern. This approach could be 
extended to include instabilities associated with modes that are tangent to the unit 
circle by looking at coupled Newell-Whitehead-Segel equations, along the lines sug- 
gested by [27], although high-order nonlinear interaction may not be captured in a 
long-wave analysis truncated at cubic order. 
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Appendix A. Weakly nonlinear theory. 

In tliis appendix, wc present tlic weakly nonlinear theory for the PDE (jS.ip . 
We will describe the calculation in terms of a harmonic primary bifurcation; the 
subharmonic case is similar, with the main differences being that the period is Att 
rather then 27r, and that the quadratic coefficient Q is identically zero. 

We start by writing U = u + iv, where u{x, y, t) and v{x, y, t) are real functions, 
and so 



ot ^ ' ^ ' 

We define differential operators £ and M: 

£ 4- - (m + aV^ + 7V*) and = (t^ + /^V^ + ,5V^) , 

at ^ ' ^ ' 

so the PDEs for u and v are 

£u = -Ali; + NL„, 
Lv = Mu + NL„ + 
The nonlinear terms NL„ and NL^, are: 

NL„ = Qiriu^ - v^) - Qi,(2ut;) + Q2r{u^ + w^) + Cr[u^ + r?-)u ~ Ci[u^ + r?-)v, 

A.l. Linear theory. We seek solutions of the form u ~ e''"'^pi(i) and v — 
e^^'^ q\(i) , where pi and qx are periodic functions of period T , and define 



du 
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we get 

Cipi = -Miqi, 

Ciqi= Mipi + f{t)pi, 

or 

Cjpi = -MiCiqi = -M\pi - f(t)MiPi, 
where Ci and A4i act on pi{t) and qi{t) as 

£i = 4- + — and Mi ^ Cli. 
dt 2 

The Hncariscd PDE reduces to a damped Mathieu equation for pi : 

(^ + y) Pi + nipi + f{t)nipi^o 

or 

Pi + JiPi + {^l + ^if{t)) Pi = = Lpi. 
The adjoint equation is: 

Pi - liPi + {^1 + ^i.f{t)) Pi = = Lpi, 
with respect to an inner product 

{9,h) = ^J^ 9{t)h{t)dt, 

with T = 2tt (harmonic case) or T = 47r (subharmonic case), so {g, hh) = (Lg, h) 

For a given value of k, seeking periodic solutions of Lp = yields an eigenvalue 
problem whose eigenvalue is the amplitude of the forcing function f{t). We use 
the method of [4] to solve this eigenvalue problem with multi-frequency forcing f{t), 
providing the critical forcing amplitude. Minimising this critical forcing amplitude 
over k yields the critical wavenumber /cc, critical forcing function fc{t), and critical 
eigenfunction pi(t). The corresponding qi{t) is determined by solving Ciqi = fipi + 

/c(t)pi. 

A. 2. Rhombs. We consider / close to fc, writing f{t) = fc{t){l + £'^^2), and 
seek small-amplitude rhombic solutions associated with two wavevectors fci and fc2 at 
the critical wavenumber: ki = k2 = kc, separated by an angle 9. We formally expand 
the solution as 

u = eui + e^U2 + e'^Ms + • • • 
= e (2i(T2)e^'=-- + 22(^2)6^^-- + c.c.) pi(i) 

+ (z?e2»^-- + z^e^^^-- -t- c.c.) P2(t) + (|zip + 1^2^) PsW 

+ (ziZ2e^('=i+'^^) " + c.c.) P4(i) + (ziZ2e^('=^"'=^) " + c.c.) psW) + O(e'), 
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with a similar expression for v in terms of Vi, V2 and 1)3, and qi, (75, where T2 
is a slow time, varying on a scale of e~^, and the fmrctions P2{i), 95 (i) are to be 
determined. The form of this solution is chosen by knowing in advance the structure 
of the nonlinear terms and the modes to be generated by them. 

Substituting these expressions for u and v into the PDE and ordering in powers 
of e, we recover, at leading order in e, the linear theory. At second order in e, we split 
the PDE into terms that go as e^^'^^ "', e^'^^^ a^^ terms without spatial dependence, and 
terms that go as e'(''i+'=2)'^ and g>-{ki-k:i)-x ^ 

The terms like e^*'"^'^ lead to equations for p2 and 

£2^2 = -X2g2+NLg, 

C2q2 = M2P2+ NLg^ + /c(t)p2, 

where the linear operators are: 

M2 = {lj- 4l3kl + 166kt) = {12, 



(2) (2) 

and the nonlinear terms NLp2 and NL^2 

Nh^^ = QiM - qf) - Qu{2pm) + Q2r{p\ + g?), 

NLg^ = Qu[pI - ql) + Qlr{2pm) + Q2^{pl + qf). 

The function (72 is eliminated by operating with C2, resulting in a second-order non- 
constant coefficient inhomogcncous linear ODE for p2 : 

{Cl +Ml+ M2.f{t)) P2 = C2'ni}^^ - X2NLg. 



P2 + 72P2 + [^l + f^2/c(i)) P2={j^ + y) ^^p2 - fi2NLg^ 

This can be solved numerically for p2 using Fourier transform methods, and q2 can 
then be found. Terms that go as e^**"^'^ result in the same equation. 

Terms without spatial dependence, and terms with spatial dependence e'('"i+'"2) '^ 
and e*('"i~'"2) a;^ result in similar equations forp3, p4 andps, but with linear operators: 

'^' = S"" = I + I' M, = ^ = n„ 

Mi=(^~ 4cos2 (^0 I3kl + 16cos'^ (^0 Sk^^ = Cl^, 



Ms^ (^uj- 4sin2 pk^ + IGsin'^ 



Skf. = n 



5, 
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and nonlinear terms: 



NLg' =NLg' =Nl2' =2NLf2'. 



Note that 74, 75, Cl^ and Cl^ depend on 6, the angle between the chosen wavevectors, 
whereas 71, 72, 73, i^i, $^2 and Cl^ do not. 

At third order in e, the problem has the following structure: 

01 2 



The nonlinear terms NL^^^ and are: 

NL^^^ = 2Qir{uiU2 - V1V2) - Qli{2uiV2 + 2u2Vi) + 2Q2r{uiU2 + V1V2) 

+ Cr{u\ + v\)ui - Ci{u\ + vl)vi, 
NL[^^ = 2Qu{uiU2 - W1W2) + (3lr(2uiW2 + 2u2Vi) + 2Q2i{uiU2 + V1V2) 

+ Cr{ul + vl)vi + C^{ul + vl)ui. 

Eliminating U3, we obtain 

(£2 + X2 ^ X/,) U3 = + - Xi^2./cui + /:NL[f) - XNL(3). 

CJ 2 01 2 

The operator on the left is the singular operator from the linearised problem, so the 
equation can only be solved for M3 if a solvability condition is applied to the terms 
that are proportional to e^^'^'^ and e**"^'^ and complex conjugates. If we take the 
inner product between pi and the e'^^^'^ component of the above, we find 

= aF2Z, + (i|zi|2 + (^^^^p + BresW)k2p) z^, 

with a similar equation for Z2, where 



T = ( pi,2 ( pi + ypi ) ) and cr = (pi, -fii/cPi 



(using Lipi = 
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A 



indcp 



Pi, + y) (2Qi,.(pip2 - gi92 - 9193) 

- 2Qu{piq2 +Piq3 + qiP2 + qiPs) 

+ 2(32r(PlP2 + P1P3 + gi92 + 9193) 
+ SCripl + ql)pi - SQipl + ql)qi) 

- Cli (2Qir{piq2 + Piqs + qiV2 + qiPs) 

+ '2Qu{piP2 +P1P3 - qiq-i - qiqs) 
+ '2Q2i{piP2 + P1P3 + 9i<?2 + qiqs) 
+ 'scM + ql)qi + ic,{p\ + ), 



71 



dt 



2Qir{piP3 - 9193) - ^Quipiqs + qips) 



+ 2(32r(piP3 + 9193) + 6C,.(pi + ql)pi - 6Ci{pl + ql)qi 

- hi (2Qir{piq3 + qiP3) + 2Qh(piP3 - qiqs) + '2Q2%{piP3 + qiqs) 

+ 6cM + ql)qi + ^c,{p\ + ), 
Pi, + y) (2Qir(piP4 +J51P5 - qiqi - qiqs) 

- 2Qu{piq4, +Piq5 + qiPi + qiPs) 

+ 2(32r(piP4 + P1P5 + qiqi + qiqb) 

- Cli(2Qir{piq4 + Piqb + qiPi + qiP5) 

+ 2QuipiP4 + P1P5 - qiqi - qiq^) 

+ 2Q2i{piPi + P1P5 + qiqi + qiqb 



For convenience, we have separated the parts of the cross-couphng coefficient that do 
not depend on the angle between the modes (-Sji^dep) those that do (i?i-es(^))- 
We discuss below how these coefficients are then scaled. 



A. 3. Hexagons. As with rhombs, we look for / close to /c, writing f{t) = 
fc{t){l + e^i^2), but now we chose three wavcvectors, fci, fc2 and ks oriented at 120° 
to each other, with fci + A;2 + fca =0. We look for small- amplitude solutions with 
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equal amplitudes of the three waves, and write 

u ~ eui + e^U2 + e'^ua + • • • 

= ez(Ti, T2) (e''^!-^ + e*'==-^ + e^'"^'^ + c.c.) pi{t) 

+ ^2 (e''=i-'^ + 6^'== =^ + e'*=^-^ + c.c.) P6(t)) + C(e^), 

with a similar expression for v in terms of gi, gg, where Ti and T2 are slow times, 
varying on scales and e~2^ and the functions P2(i), Qeit) a-re to be determined. 
The form of the expression is chosen by knowing in advance the structure of the 
nonlinear terms. 

Substituting these expressions for u and v into the PDE and ordering in powers 
of e, at leading order in e we recover the linear theory. At second order in e, we split 
the PDE into terms that go as e2''"i'^ + e2''^2 a; _|_ g2i/c3 x _j_ ^ terms without spatial 
dependence, terms that go as e*('=i^'=2)-^ + Q^ik2-ks} x _|_ gi{fc3-fci)-x _|_ ^ finally 
terms that go as e'^^i'^ + e^'^^'^ + e^'^^'^ + c.c, which have to be considered specially. 

Terms like e^^^'^'^ + e2*'"2'^ + e"^^^^'^ + c.c. and terms without spatial dependence 
give the same equations for p2, 92, and as in the case of rhombs. In particular, 
the inhomogeneous nonlinear terms are the same. Terms that go as e'^'^^-'^^) ^ _^ 
gi(fc2-fc3) x _|_ ^i(k3-ki)-x result in equations for p^ and ^5 that are the equations 

for p5 and 95 evaluated for 6 = 120°. 

Terms that go as e^^'^'^ + e^^^ '^ + e^'^^'^ + c.c. require the use of two time-scales 
and a solvability condition. The linear operators are the same as for the initial linear 
problem: 

,A , ATT (2) d^'/dTl 

Lipe = -MiciQ + NL^g' '-^ — pi, 

£196= Mip^ + n\}^^ + ,f,lt)p^-^^^qi. 

(21 (2) 

The nonlinear terms are NL^g and NL^g are: 

NL^'g) = 2NL|,'2) and NL^g^ = 2NLf2^ 

This can be reduced to a second-order non-constant coefficient inhomogeneous linear 
ODE for Pa: 

{C\ +Ml + Mif{t)) pe = CNlI^^ - XiNLfg) + {Mm - Cm) 

Since the operator on the LHS is the singular linear operator L, we must apply a 
solvability condition: 
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since MiQi — —Cipi. We define 



as before and 



~ dpi 



e= (pi,/:iNLg^-AliNLfe) 



and obtain an equation for the slow (Ti) evolution of the amplitude z: 

dz 



Once the solvability condition has been imposed, the ODE can be solved for pg and 
Qq. The computed solution pg contains an arbitrary amount of pi; the solution is 
made unique by specifying that (pi^pe) = 0. 

At third order in e, the problem has the following structure: 

oil OI2 

We only need keep track of terms proportional to e'^^^ '^ in our derivation of the bifur- 
cation problem, so the duijdTi and dv^jdTi terms yield p^dz^ / dTi and qedz'^ /dTi. 
The e*'^! '^ components of the nonlinear terms and are specified below. 

Eliminating V3, we obtain 

(£- + A1» + Mf.) „3 = -£ + ll - NL<»»: 

The equation can only be solved for U3 if a solvability condition is applied to the 
terms that are proportional to e**'^'^, e^^^'^ and e**"^'^, and complex conjugates. If 
we take the inner product between pi and the e*'"^'^ component of the above, we find 

- aF2Z +{a + 2Bgo) \z\^z, 



dT2 

where r and a are unchanged from the rhombic calculations, and 
A + 2B,o = (pi, (I + |)nl(3'') - fiiNLl^") 
+ (-Cipe+Miqe)2^y 

where we have used 
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and 

NL^^''^ = 2Qir{piP2 + 3piP3 + 2pip5 + 2pipe - qiq2 - SqiQs - 217195 - 2qiqe) 

- 2Qii{piq2 + 3piq3 + 2piq5 + 2piqe + qiP2 + 'SqiP3 + 2qipz + 2qipQ) 
+ 2Qir{piP2 + 3piP3 + 2pipz + 2pipe + qiq2 + 'iqiqa + 2(?iQ5 + 2qiqe) 
+ IdCripj + ql)Pi - 15a(p? + qi)qu 

NL^^''^ = 2Qir{piq2 + 3piq3 + 2piq5 + 2piqe + qiP2 + iqipa + 2qip5 + 2qipe) 

+ 2Qu{pip2 + ipiPa + 2piP5 + 2pip6 - qiq2 - iqiqs - 217195 - 2qiqe) 
+ 2Q2i{piP2 + ipiPa + 2piP5 + 2pipe + qiq2 + iqiqa + 2qiq^ + 2qiqe) 

+ IbCripl + qDQI + Cripl + qDpI- 

From the value of A calculated for rhombs, we can recover BgOi which is effectively 
the cross-coupling coefficient for modes with wavevectors at 60°. 

A. 4. Reconstitution. At this stage, the pattern formation problem on a hexag- 
onal lattice would take the form: 
dzi 

= '''''' 

= <7F2Zi + (i|zip + Beo\z2\^ + Beolzsl"") zi, 

where Zi, Z2 and Z3 are amplitudes of e^'^^'^, Qikz-x ^iks-x gjjjjjjj^j. equations 
are found for dz2/dTi etc. Recall that the small factor e has been used so that the 
amplitude of the original amplitude U is explicitly small: U — €Zx{Tx,T2)e'^'^^'^p\{t) + 

There is more than one way to combine these equations into a single ODE. Prop- 
erly, we should consider only the case where the coefficient of the quadratic term e 
is itself small (order e). This occurs either for small values of Qi and Q2, or near a 
codimension-one line in {Qi, Q2) space: 

= e = (pi, /IiNL^'g) - A^iNLfg)) 

= 2Qi, (piXiipl - qi) - Mi{2pm)) + 2Qu {pi,£i{--2pm) - MM - g?)) 
+ 2Q2r {pi,Ci{pl + ql))+ 2Q2. {pi,~MM + qi)). 

Alternatively, we note that for many of these multi-frequency forced problems, e is 
small anyway [47,48]. 

Having decided that e is small, we define a fast time scale t* , related to the original 
time scale t by an order-one factor cr/r: 

d T d 2'^ 
lit* ^ ^'^Wi^^ adf^' 
where e is the original small parameter. We now scale the z's by 1/e. so that the 
original amplitude U is implicitly small: U = zi{Ti,T2)e^'^^''^pi{t) -I- . . . , and obtain: 

^^Xzi + -Z2Z3 + - (|2;i|2 + Sgol^aP + Sgolzap) zi, 
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where A = e^i^2i so that the forcing amphtude is (1 + A) times the critical amphtude, 



The advantage of reconstituting in this way is that the quadratic and cubic terms 
appear at the same order. The disadvantage is that the regime of vahdity (A ^ 1, 
e ^ 1, 2; ^ 1) is not made exphcit. In particular, this validity condition will only be 
satisfied for hexagons when the coefficient of the quadratic term is small - which is 
precisely the limit required for the quadratic and cubic terms to be of the same order. 

Finally, we scale the amplitudes once more by a factor of \J \(J / A\, rename the 
time variable back to t, and obtain: 



and Bi 



'60 



B^a/A. 



dzi 
~dt 



Azi + Qz2Z3 + s (|zip + Bqo\z2\^ + BeokaP) zi 



with similar equations for Z2 and Z3, where 




and 



s = sgn 




(usually, s = —1). 

Repeating the same reconstitution for the rhombic lattice results in 



dzi 
~dt 
dz2 
~dt 



Azi +s(|zip + B(,|z2p)zi 



AZ2+S(|Z2P+B(,|Z1P)Z2 



where Bg = (Bjj^j^p + A-es(f))/i- 

The equations for other lattices can be found from combinations of the above, as 
can candidate equations for quasipatterns. 



